gwf_solver.F90

Go to the documentation of this file.

!     ###############################################################################
      SUBROUTINE gwf_solver (KLON,KLAT,PNPTS,OMASK,PHCOF,PRHS,PCR,PCC,PHGROUND,PEVOL)
!     ###############################################################################
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
INTEGER, INTENT(IN)                 :: KLON
INTEGER, INTENT(IN)                 :: KLAT
!
REAL, INTENT(IN)                    :: PNPTS
!
LOGICAL, DIMENSION(:,:), INTENT(IN) :: OMASK
!
REAL, DIMENSION(:,:), INTENT(IN)    :: PHCOF
REAL, DIMENSION(:,:), INTENT(IN)    :: PRHS
REAL, DIMENSION(:,:), INTENT(IN)    :: PCR
REAL, DIMENSION(:,:), INTENT(IN)    :: PCC
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PHGROUND
REAL,                 INTENT(OUT)   :: PEVOL
!
!*      0.2    declarations of local variables
!
INTEGER                        :: JLON, JLAT
INTEGER                        :: INBSOU
REAL                           :: ZHOR,ZNOM,ZDENOM,ZDIAG,ZPASD,ZDIAG_HORI
REAL                           :: ZEVOL
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
! *      1. INITIALISATION
!           --------------
!
IF (lhook) CALL dr_hook('GWF_SOLVER',0,zhook_handle)
!
! * error initialization
!
zevol = 0.
!
! *       2. ONE ITERATION LOOP
!            ------------------
!
DO jlat = 1,klat
    DO jlon = 1,klon
!
        IF (omask(jlon,jlat)) THEN
!
!          skip constant head cells and no-flow
           zhor   = 0.
           zdiag  = 0.
           znom   = 0.
           zdenom = 0.
           zpasd  = 0.
!
!          RIGHT FLUX
           IF (jlon < klon) THEN
              zdiag_hori = pcr(jlon,jlat)
              zhor  = zhor  + zdiag_hori*phground(jlon+1,jlat)
              zdiag = zdiag + zdiag_hori
           ENDIF
!       
!          LEFT FLUX
           IF (jlon > 1) THEN
              zdiag_hori = pcr(jlon-1,jlat)
              zhor  = zhor  + zdiag_hori*phground(jlon-1,jlat)
              zdiag = zdiag + zdiag_hori
           ENDIF
!           
!          BOTTOM FLUX
           IF (jlat > 1) THEN
              zdiag_hori = pcc(jlon,jlat-1)
              zhor  = zhor  + zdiag_hori*phground(jlon,jlat-1)
              zdiag = zdiag + zdiag_hori
           ENDIF
!           
!          TOP FLUX
           IF (jlat < klat) THEN
              zdiag_hori = pcc(jlon,jlat)
              zhor  = zhor  + zdiag_hori*phground(jlon,jlat+1)
              zdiag = zdiag + zdiag_hori
           ENDIF
!
!          calcul nominateur
           znom   = zhor + prhs(jlon,jlat)
!           
!          calcul dénominateur
           zdenom = zdiag + phcof(jlon,jlat)
!
!          calcul new height
           zpasd = znom/zdenom
!        
!          calcul critère convergence
           zevol = zevol + abs(zpasd-phground(jlon,jlat))/pnpts
!
!          update new height
           phground(jlon,jlat) = zpasd
!
        ENDIF
!
    ENDDO
ENDDO
!
pevol = zevol
!
IF (lhook) CALL dr_hook('GWF_SOLVER',1,zhook_handle)
!
END SUBROUTINE gwf_solver