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 5842 for branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 – NEMO

Ignore:
Timestamp:
2015-10-30T16:34:24+01:00 (8 years ago)
Author:
hliu
Message:

Wetting and Drying update based on r:5803

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5224 r5842  
    1717   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
    1818   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity 
     19   !!            3.6? !  2015-11  (H. Liu) add Wetting/Drying pressure filter 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    378379      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    379380      !! 
    380       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    381       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     381      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
     382      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars 
     383      LOGICAL  ::   ll_tmp1, ll_tmp2, ll_tmp3            ! local logical variables 
    382384      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     385      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter 
    383386      !!---------------------------------------------------------------------- 
    384387      ! 
    385388      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     389      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    386390      ! 
    387391      IF( kt == nit000 ) THEN 
     
    397401      ELSE                 ;     znad = 0._wp         ! Fixed volume 
    398402      ENDIF 
     403 
     404      IF(ln_wd) THEN 
     405        DO jj = 2, jpjm1 
     406           DO ji = 2, jpim1  
     407             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))  
     408             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 
     409             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     410                                                       & rn_wdmin1 + rn_wdmin2 
     411 
     412             IF(ll_tmp1.AND.ll_tmp2) THEN 
     413               zcpx(ji,jj) = 1.0_wp 
     414               wduflt(ji,jj) = 1.0_wp 
     415             ELSE IF(ll_tmp3) THEN 
     416               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     417               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 
     418                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     419               wduflt(ji,jj) = 1.0_wp 
     420             ELSE 
     421               zcpx(ji,jj) = 0._wp 
     422               wduflt(ji,jj) = 0.0_wp 
     423             END IF 
     424       
     425             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))  
     426             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 
     427             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     428                                                       & rn_wdmin1 + rn_wdmin2 
     429 
     430             IF(ll_tmp1.AND.ll_tmp2) THEN 
     431               zcpy(ji,jj) = 1.0_wp 
     432               wdvflt(ji,jj) = 1.0_wp 
     433             ELSE IF(ll_tmp3) THEN 
     434               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     435               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 
     436                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     437               wdvflt(ji,jj) = 1.0_wp 
     438             ELSE 
     439               zcpy(ji,jj) = 0._wp 
     440               wdvflt(ji,jj) = 0.0_wp 
     441             END IF 
     442           END DO 
     443        END DO 
     444        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     445      ENDIF 
     446 
     447 
     448      !jii=jidbg-nimpp+1;jjj=jjdbg-njmpp+1 
    399449 
    400450      ! Surface value 
     
    411461            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    412462               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     463 
     464 
     465            IF(ln_wd) THEN 
     466 
     467              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     468              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     469              zuap = zuap * zcpx(ji,jj) 
     470              zvap = zvap * zcpy(ji,jj) 
     471            ENDIF 
     472 
    413473            ! add to the general momentum trend 
    414474            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    433493               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    434494                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     495 
     496               IF(ln_wd) THEN 
     497                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     498                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     499                 zuap = zuap * zcpx(ji,jj) 
     500                 zvap = zvap * zcpy(ji,jj) 
     501               ENDIF 
     502 
    435503               ! add to the general momentum trend 
    436504               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    441509      ! 
    442510      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     511      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    443512      ! 
    444513   END SUBROUTINE hpg_sco 
     
    719788      REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    720789      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     790      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    721791      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    722792      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw 
    723793      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow 
    724794      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k 
     795      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    725796      !!---------------------------------------------------------------------- 
    726797      ! 
     
    728799      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    729800      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
    730       ! 
     801      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     802      ! 
     803      !!---------------------------------------------------------------------- 
     804      ! 
     805      ! 
     806      IF(ln_wd) THEN 
     807        DO jj = 2, jpjm1 
     808           DO ji = 2, jpim1  
     809             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     810                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     811                     &  > rn_wdmin1 + rn_wdmin2 
     812             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     813                                                       & rn_wdmin1 + rn_wdmin2 
     814 
     815             IF(ll_tmp1) THEN 
     816               zcpx(ji,jj) = 1.0_wp 
     817             ELSE IF(ll_tmp2) THEN 
     818               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     819               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     820                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     821             ELSE 
     822               zcpx(ji,jj) = 0._wp 
     823             END IF 
     824       
     825             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     826                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     827                     &  > rn_wdmin1 + rn_wdmin2 
     828             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     829                                                       & rn_wdmin1 + rn_wdmin2 
     830 
     831             IF(ll_tmp1) THEN 
     832               zcpy(ji,jj) = 1.0_wp 
     833             ELSE IF(ll_tmp2) THEN 
     834               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     835               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     836                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     837             ELSE 
     838               zcpy(ji,jj) = 0._wp 
     839             END IF 
     840           END DO 
     841        END DO 
     842        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     843      ENDIF 
     844 
    731845 
    732846      IF( kt == nit000 ) THEN 
     
    8991013            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 
    9001014            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 
     1015            IF(ln_wd) THEN 
     1016              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     1017              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     1018            ENDIF 
    9011019            ! add to the general momentum trend 
    9021020            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    9181036                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    9191037                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj) 
     1038               IF(ln_wd) THEN 
     1039                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     1040                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     1041               ENDIF 
    9201042               ! add to the general momentum trend 
    9211043               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    9281050      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    9291051      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     1052      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    9301053      ! 
    9311054   END SUBROUTINE hpg_djc 
     
    9511074      !! The local variables for the correction term 
    9521075      INTEGER  :: jk1, jis, jid, jjs, jjd 
     1076      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables 
    9531077      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
    9541078      REAL(wp) :: zrhdt1 
     
    9571081      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
    9581082      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n 
     1083      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    9591084      !!---------------------------------------------------------------------- 
    9601085      ! 
     
    9621087      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 
    9631088      CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 
     1089      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    9641090      ! 
    9651091      IF( kt == nit000 ) THEN 
     
    9741100      znad = 0.0_wp 
    9751101      IF( lk_vvl ) znad = 1._wp 
     1102 
     1103      IF(ln_wd) THEN 
     1104        DO jj = 2, jpjm1 
     1105           DO ji = 2, jpim1  
     1106             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     1107                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     1108                     &  > rn_wdmin1 + rn_wdmin2 
     1109             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     1110                                                       & rn_wdmin1 + rn_wdmin2 
     1111 
     1112             IF(ll_tmp1) THEN 
     1113               zcpx(ji,jj) = 1.0_wp 
     1114             ELSE IF(ll_tmp2) THEN 
     1115               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     1116               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1117                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     1118             ELSE 
     1119               zcpx(ji,jj) = 0._wp 
     1120             END IF 
     1121       
     1122             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     1123                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     1124                     &  > rn_wdmin1 + rn_wdmin2 
     1125             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     1126                                                       & rn_wdmin1 + rn_wdmin2 
     1127 
     1128             IF(ll_tmp1.OR.ll_tmp2) THEN 
     1129               zcpy(ji,jj) = 1.0_wp 
     1130             ELSE IF(ll_tmp2) THEN 
     1131               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     1132               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1133                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     1134             ELSE 
     1135               zcpy(ji,jj) = 0._wp 
     1136             END IF 
     1137           END DO 
     1138        END DO 
     1139        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     1140      ENDIF 
    9761141 
    9771142      ! Clean 3-D work arrays 
     
    10521217        END DO 
    10531218      END DO 
     1219 
     1220      CALL lbc_lnk (zsshu_n, 'U', 1.) 
     1221      CALL lbc_lnk (zsshv_n, 'V', 1.) 
    10541222 
    10551223      DO jj = 2, jpjm1 
     
    11501318               ENDIF 
    11511319 
    1152                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 
    1153                &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1320               IF(ln_wd) THEN 
     1321                  zdpdx1 = zdpdx1 * zcpx(ji,jj) 
     1322                  zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1323                ENDIF 
     1324                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    11541325            ENDIF 
    11551326 
     
    12071378               ENDIF 
    12081379 
    1209                va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 
    1210                &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1380               IF(ln_wd) THEN 
     1381                  zdpdy1 = zdpdy1 * zcpy(ji,jj) 
     1382                  zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1383                ENDIF 
     1384 
     1385               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    12111386            ENDIF 
    12121387 
     
    12191394      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 
    12201395      CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 
     1396      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    12211397      ! 
    12221398   END SUBROUTINE hpg_prj 
Note: See TracChangeset for help on using the changeset viewer.