New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6152 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 – NEMO

Ignore:
Timestamp:
2015-12-21T23:33:57+01:00 (8 years ago)
Author:
acc
Message:

Add wetting and drying option from dev_r5803_NOC_WAD branch. Logically isolated code changes in domvvl.F90, domzgr.F90, dynhpg.F90, dynspg_ts.F90, sshwzv.F90 and nemogcm.F90. New module wet_dry.F90 in DYN. Fully SETTE tested with code deactivated (ln_wad=.false.). No test case yet available to justify activating option (still under development)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r6140 r6152  
    3333   USE sbc_oce         ! surface variable (only for the flag with ice shelf) 
    3434   USE dom_oce         ! ocean space and time domain 
     35   USE wet_dry         ! wetting and drying 
    3536   USE phycst          ! physical constants 
    3637   USE trd_oce         ! trends: ocean variables 
    … …  
    429430      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    430431      !! 
    431       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    432       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     432      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
     433      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars 
     434      LOGICAL  ::   ll_tmp1, ll_tmp2, ll_tmp3            ! local logical variables 
    433435      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     436      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter 
    434437      !!---------------------------------------------------------------------- 
    435438      ! 
    436439      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     440      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    437441      ! 
    438442      IF( kt == nit000 ) THEN 
    … …  
    447451      ENDIF 
    448452      ! 
     453      IF(ln_wd) THEN 
     454        DO jj = 2, jpjm1 
     455           DO ji = 2, jpim1  
     456             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))  
     457             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 
     458             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     459                                                       & rn_wdmin1 + rn_wdmin2 
     460 
     461             IF(ll_tmp1.AND.ll_tmp2) THEN 
     462               zcpx(ji,jj) = 1.0_wp 
     463               wduflt(ji,jj) = 1.0_wp 
     464             ELSE IF(ll_tmp3) THEN 
     465               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     466               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 
     467                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     468               wduflt(ji,jj) = 1.0_wp 
     469             ELSE 
     470               zcpx(ji,jj) = 0._wp 
     471               wduflt(ji,jj) = 0.0_wp 
     472             END IF 
     473       
     474             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))  
     475             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 
     476             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     477                                                       & rn_wdmin1 + rn_wdmin2 
     478 
     479             IF(ll_tmp1.AND.ll_tmp2) THEN 
     480               zcpy(ji,jj) = 1.0_wp 
     481               wdvflt(ji,jj) = 1.0_wp 
     482             ELSE IF(ll_tmp3) THEN 
     483               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     484               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 
     485                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     486               wdvflt(ji,jj) = 1.0_wp 
     487             ELSE 
     488               zcpy(ji,jj) = 0._wp 
     489               wdvflt(ji,jj) = 0.0_wp 
     490             END IF 
     491           END DO 
     492        END DO 
     493        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     494      ENDIF 
     495 
     496 
    449497      ! Surface value 
    450498      DO jj = 2, jpjm1 
    … …  
    460508            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    461509               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
     510 
     511 
     512            IF(ln_wd) THEN 
     513 
     514              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     515              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     516              zuap = zuap * zcpx(ji,jj) 
     517              zvap = zvap * zcpy(ji,jj) 
     518            ENDIF 
     519 
    462520            ! add to the general momentum trend 
    463521            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    … …  
    482540               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
    483541                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 
     542 
     543               IF(ln_wd) THEN 
     544                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     545                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     546                 zuap = zuap * zcpx(ji,jj) 
     547                 zvap = zvap * zcpy(ji,jj) 
     548               ENDIF 
     549 
    484550               ! add to the general momentum trend 
    485551               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    … …  
    490556      ! 
    491557      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     558      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    492559      ! 
    493560   END SUBROUTINE hpg_sco 
    … …  
    623690      REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    624691      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     692      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    625693      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    626694      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw 
    627695      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow 
    628696      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k 
     697      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    629698      !!---------------------------------------------------------------------- 
    630699      ! 
    … …  
    632701      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    633702      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
    634       ! 
     703      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     704      ! 
     705      ! 
     706      IF(ln_wd) THEN 
     707        DO jj = 2, jpjm1 
     708           DO ji = 2, jpim1  
     709             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     710                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     711                     &  > rn_wdmin1 + rn_wdmin2 
     712             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     713                                                       & rn_wdmin1 + rn_wdmin2 
     714 
     715             IF(ll_tmp1) THEN 
     716               zcpx(ji,jj) = 1.0_wp 
     717             ELSE IF(ll_tmp2) THEN 
     718               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     719               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     720                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     721             ELSE 
     722               zcpx(ji,jj) = 0._wp 
     723             END IF 
     724       
     725             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     726                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     727                     &  > rn_wdmin1 + rn_wdmin2 
     728             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     729                                                       & rn_wdmin1 + rn_wdmin2 
     730 
     731             IF(ll_tmp1) THEN 
     732               zcpy(ji,jj) = 1.0_wp 
     733             ELSE IF(ll_tmp2) THEN 
     734               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     735               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     736                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     737             ELSE 
     738               zcpy(ji,jj) = 0._wp 
     739             END IF 
     740           END DO 
     741        END DO 
     742        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     743      ENDIF 
     744 
    635745 
    636746      IF( kt == nit000 ) THEN 
    … …  
    803913            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
    804914            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
     915            IF(ln_wd) THEN 
     916              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     917              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     918            ENDIF 
    805919            ! add to the general momentum trend 
    806920            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
    … …  
    822936                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    823937                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
     938               IF(ln_wd) THEN 
     939                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     940                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     941               ENDIF 
    824942               ! add to the general momentum trend 
    825943               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    … …  
    832950      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    833951      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     952      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    834953      ! 
    835954   END SUBROUTINE hpg_djc 
    … …  
    855974      !! The local variables for the correction term 
    856975      INTEGER  :: jk1, jis, jid, jjs, jjd 
     976      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables 
    857977      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
    858978      REAL(wp) :: zrhdt1 
    … …  
    861981      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
    862982      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n 
     983      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    863984      !!---------------------------------------------------------------------- 
    864985      ! 
    … …  
    866987      CALL wrk_alloc( jpi,jpj,jpk,   zdept, zrhh ) 
    867988      CALL wrk_alloc( jpi,jpj,       zsshu_n, zsshv_n ) 
     989      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    868990      ! 
    869991      IF( kt == nit000 ) THEN 
    … …  
    877999      znad = 1._wp 
    8781000      IF( ln_linssh )   znad = 0._wp 
     1001 
     1002      IF(ln_wd) THEN 
     1003        DO jj = 2, jpjm1 
     1004           DO ji = 2, jpim1  
     1005             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     1006                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     1007                     &  > rn_wdmin1 + rn_wdmin2 
     1008             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     1009                                                       & rn_wdmin1 + rn_wdmin2 
     1010 
     1011             IF(ll_tmp1) THEN 
     1012               zcpx(ji,jj) = 1.0_wp 
     1013             ELSE IF(ll_tmp2) THEN 
     1014               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     1015               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1016                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     1017             ELSE 
     1018               zcpx(ji,jj) = 0._wp 
     1019             END IF 
     1020       
     1021             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     1022                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     1023                     &  > rn_wdmin1 + rn_wdmin2 
     1024             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     1025                                                       & rn_wdmin1 + rn_wdmin2 
     1026 
     1027             IF(ll_tmp1.OR.ll_tmp2) THEN 
     1028               zcpy(ji,jj) = 1.0_wp 
     1029             ELSE IF(ll_tmp2) THEN 
     1030               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     1031               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1032                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     1033             ELSE 
     1034               zcpy(ji,jj) = 0._wp 
     1035             END IF 
     1036           END DO 
     1037        END DO 
     1038        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     1039      ENDIF 
    8791040 
    8801041      ! Clean 3-D work arrays 
    … …  
    9601121        END DO 
    9611122      END DO 
     1123 
     1124      CALL lbc_lnk (zsshu_n, 'U', 1.) 
     1125      CALL lbc_lnk (zsshv_n, 'V', 1.) 
    9621126 
    9631127      DO jj = 2, jpjm1 
    … …  
    10571221                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
    10581222               ENDIF 
    1059 !!gm  Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:)  by definition 
    1060 !!gm      in the line below only * umask(ji,jj,jk)  is needed !! 
    1061                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1223               IF(ln_wd) THEN 
     1224                  zdpdx1 = zdpdx1 * zcpx(ji,jj) 
     1225                  zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1226               ENDIF 
     1227               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    10621228            ENDIF 
    10631229 
    … …  
    11141280                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
    11151281               ENDIF 
    1116 !!gm  Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:)  by definition 
    1117 !!gm      in the line below only * vmask(ji,jj,jk)  is needed !! 
    1118                va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1282               IF(ln_wd) THEN 
     1283                  zdpdy1 = zdpdy1 * zcpy(ji,jj) 
     1284                  zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1285               ENDIF 
     1286 
     1287               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    11191288            ENDIF 
    11201289               ! 
    … …  
    11261295      CALL wrk_dealloc( jpi,jpj,jpk,   zdept, zrhh ) 
    11271296      CALL wrk_dealloc( jpi,jpj,       zsshu_n, zsshv_n ) 
     1297      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    11281298      ! 
    11291299   END SUBROUTINE hpg_prj 
Note: See TracChangeset for help on using the changeset viewer.