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 9023 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 – NEMO

Ignore:
Timestamp:
2017-12-13T18:08:50+01:00 (6 years ago)
Author:
timgraham
Message:

Merged METO_MERCATOR branch and resolved all conflicts in OPA_SRC

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90

    r9019 r9023  
    11MODULE wet_dry 
     2 
     3   !! includes updates to namelist namwad for diagnostic outputs of ROMS wetting and drying 
     4 
    25   !!============================================================================== 
    36   !!                       ***  MODULE  wet_dry  *** 
    47   !! Wetting and drying includes initialisation routine and routines to 
    58   !! compute and apply flux limiters and preserve water depth positivity 
    6    !! only effects if wetting/drying is on (ln_wd == .true.) 
     9   !! only effects if wetting/drying is on (ln_wd_il == .true. or ln_wd_dl==.true. ) 
    710   !!============================================================================== 
    811   !! History :  3.6  ! 2014-09  ((H.Liu)  Original code 
     
    3336 
    3437   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask   !: u- and v- limiter  
    35    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd    !: wetting and drying t-pt depths 
    3638   !                                                           !  (can include negative depths) 
    37  
    38    LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off) 
     39   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdramp, wdrampu, wdrampv !: for hpg limiting 
     40 
     41   LOGICAL,  PUBLIC  ::   ln_wd_il    !: Wetting/drying il activation switch (T:on,F:off) 
     42   LOGICAL,  PUBLIC  ::   ln_wd_dl    !: Wetting/drying dl activation switch (T:on,F:off) 
     43   REAL(wp), PUBLIC  ::   rn_wdmin0   !: depth at which wetting/drying starts 
    3944   REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells 
    40    REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells 
     45   REAL(wp), PUBLIC  ::   r_rn_wdmin1 !: 1/minimum water depth on dried cells  
     46   REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerance of minimum water depth on dried cells 
    4147   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying will be considered 
    4248   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter 
     49   LOGICAL,  PUBLIC  ::   ln_wd_dl_bc !: DL scheme: True implies 3D velocities are set to the barotropic values at points  
     50                                      !: where the flow is from wet points on less than half the barotropic sub-steps   
     51   LOGICAL,  PUBLIC  ::  ln_wd_dl_rmp !: use a ramp for the dl flux limiter between 2 rn_wdmin1 and rn_wdmin1 (rather than a cut-off at rn_wdmin1)       
     52   REAL(wp), PUBLIC  ::   ssh_ref     !: height of z=0 with respect to the geoid;  
     53 
     54   LOGICAL,  PUBLIC  ::   ll_wd       !: Wetting/drying activation switch if either ln_wd_il or ln_wd_dl 
    4355 
    4456   PUBLIC   wad_init                  ! initialisation routine called by step.F90 
     
    5971      !! ** input   : - namwad namelist 
    6072      !!---------------------------------------------------------------------- 
    61       INTEGER  ::   ios, ierr   ! Local integer 
    62       !! 
    63       NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 
     73      !! 
     74      NAMELIST/namwad/ ln_wd_il, ln_wd_dl, rn_wdmin0, rn_wdmin1, rn_wdmin2, rn_wdld, & 
     75                     & nn_wdit, ln_wd_dl_bc, ln_wd_dl_rmp 
     76      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     77      INTEGER  ::   ierr                ! Local integer status array allocation  
    6478      !!---------------------------------------------------------------------- 
    6579      ! 
     
    7791         WRITE(numout,*) '~~~~~~~~' 
    7892         WRITE(numout,*) '   Namelist namwad' 
    79          WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd 
     93         WRITE(numout,*) '      Logical for Iter Lim wd option   ln_wd_il     = ', ln_wd_il 
     94         WRITE(numout,*) '      Logical for Dir. Lim wd option   ln_wd_dl     = ', ln_wd_dl 
     95         WRITE(numout,*) '      Depth at which wet/drying starts rn_wdmin0    = ', rn_wdmin0 
    8096         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1 
    81          WRITE(numout,*) '      Tolerance of min wet depth     rn_wdmin2      = ', rn_wdmin2 
     97         WRITE(numout,*) '      Tolerance of min wet depth       rn_wdmin2    = ', rn_wdmin2 
    8298         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld 
    8399         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit 
     100         WRITE(numout,*) '      T => baroclinic u,v=0 at dry pts: ln_wd_dl_bc = ', ln_wd_dl_bc      
     101         WRITE(numout,*) '      use a ramp for rwd limiter:  ln_wd_dl_rwd_rmp = ', ln_wd_dl_rmp 
     102  
    84103      ENDIF 
    85       ! 
    86       IF(ln_wd) THEN 
    87          ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj),  STAT=ierr ) 
     104      IF( .NOT. ln_read_cfg ) THEN 
     105         WRITE(numout,*) '      No configuration file so seting ssh_ref to zero  ' 
     106          ssh_ref=0.0 
     107      ENDIF 
     108 
     109      r_rn_wdmin1=1/rn_wdmin1 
     110      ll_wd = .FALSE. 
     111      IF(ln_wd_il .OR. ln_wd_dl) THEN 
     112         ll_wd = .TRUE. 
     113         ALLOCATE( wdmask(jpi,jpj),   STAT=ierr ) 
     114         ALLOCATE( wdramp(jpi,jpj), wdrampu(jpi,jpj), wdrampv(jpi,jpj), STAT=ierr )  
    88115         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 
    89116      ENDIF 
     
    118145      !!---------------------------------------------------------------------- 
    119146      ! 
    120       IF( ln_timing )   CALL timing_start('wad_lmt') 
    121       ! 
    122       !IF(lwp) WRITE(numout,*) 
    123       !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting' 
    124       ! 
     147      IF( nn_timing == 1 )  CALL timing_start('wad_lmt') 
     148      ! 
     149        
     150      DO jk = 1, jpkm1 
     151         un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:)  
     152         vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:,:)  
     153      END DO 
    125154      jflag  = 0 
    126155      zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen 
    127       !  
     156 
    128157      zflxp(:,:)   = 0._wp 
    129158      zflxn(:,:)   = 0._wp 
    130159      zflxu(:,:)   = 0._wp 
    131160      zflxv(:,:)   = 0._wp 
    132       ! 
    133       zwdlmtu(:,:) = 1._wp 
    134       zwdlmtv(:,:) = 1._wp 
    135       !  
     161 
     162      zwdlmtu(:,:)  = 1._wp 
     163      zwdlmtv(:,:)  = 1._wp 
     164       
    136165      ! Horizontal Flux in u and v direction 
    137166      DO jk = 1, jpkm1   
     
    143172         END DO   
    144173      END DO 
    145       ! 
     174        
    146175      zflxu(:,:) = zflxu(:,:) * e2u(:,:) 
    147176      zflxv(:,:) = zflxv(:,:) * e1v(:,:) 
    148       !  
     177       
    149178      wdmask(:,:) = 1 
    150179      DO jj = 2, jpj 
    151180         DO ji = 2, jpi  
    152             ! 
    153             IF( tmask(ji, jj, 1) < 0.5_wp )   CYCLE   ! we don't care about land cells 
    154             IF( ht_wd(ji,jj)     > zdepwd )  CYCLE   ! and cells which are unlikely to dry 
    155             ! 
    156             zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj) , 0._wp )  & 
    157                &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,jj-1) , 0._wp )  
    158             zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj) , 0._wp )  & 
    159                &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,jj-1) , 0._wp )  
    160             ! 
    161             zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 
    162             IF( zdep2 .le. 0._wp) THEN  !add more safty, but not necessary 
    163                sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
    164                IF( zflxu(ji,  jj) > 0._wp )  zwdlmtu(ji  ,jj) = 0._wp 
    165                IF( zflxu(ji-1,jj) < 0._wp )  zwdlmtu(ji-1,jj) = 0._wp 
    166                IF( zflxv(ji,  jj) > 0._wp )  zwdlmtv(ji  ,jj) = 0._wp 
    167                IF( zflxv(ji,jj-1) < 0._wp )  zwdlmtv(ji,jj-1) = 0._wp  
     181 
     182            IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE              ! we don't care about land cells 
     183            IF( ht_0(ji,jj) - ssh_ref > zdepwd ) CYCLE   ! and cells which are unlikely to dry 
     184 
     185            zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + & 
     186                         & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp)  
     187            zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + & 
     188                         & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp)  
     189 
     190            zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 
     191            IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary 
     192               sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 
     193               IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp 
     194               IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 
     195               IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp 
     196               IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp  
    168197               wdmask(ji,jj) = 0._wp 
    169             ENDIF 
     198            END IF 
     199         ENDDO 
     200      END DO 
     201 
     202 
     203! HPG limiter from jholt 
     204      wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp) 
     205!jth assume don't need a lbc_lnk here 
     206      DO jj = 1, jpjm1 
     207         DO ji = 1, jpim1 
     208            wdrampu(ji,jj) = min(wdramp(ji,jj),wdramp(ji+1,jj)) 
     209            wdrampv(ji,jj) = min(wdramp(ji,jj),wdramp(ji,jj+1)) 
    170210         END DO 
    171211      END DO 
    172       !! 
    173       !! start limiter iterations  
     212! end HPG limiter 
     213 
     214 
     215       
     216        !! start limiter iterations  
    174217      DO jk1 = 1, nn_wdit + 1 
    175          !  
     218        
     219           
    176220         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 
    177221         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    178222         jflag = 0     ! flag indicating if any further iterations are needed 
    179          !  
     223           
    180224         DO jj = 2, jpj 
    181225            DO ji = 2, jpi  
    182                ! 
    183                IF( tmask(ji,jj,1) < 0.5_wp )  CYCLE  
    184                IF( ht_wd(ji,jj)   > zdepwd )   CYCLE 
    185                ! 
     226         
     227               IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE  
     228               IF( ht_0(ji,jj) > zdepwd )      CYCLE 
     229         
    186230               ztmp = e1e2t(ji,jj) 
    187                ! 
    188                zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj) , 0._wp )  & 
    189                   &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,jj-1) , 0._wp )  
    190                zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj) , 0._wp )  & 
    191                   &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,jj-1) , 0._wp )  
    192                ! 
     231 
     232               zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + & 
     233                      & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp)  
     234               zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + & 
     235                      & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp)  
     236           
    193237               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
    194                zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 
    195                ! 
     238               zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 
     239           
    196240               IF( zdep1 > zdep2 ) THEN 
    197241                  wdmask(ji, jj) = 0 
     
    201245                  ! changes have zeroed the coefficient since further iterations will 
    202246                  ! not change anything 
    203                   IF( zcoef > 0._wp ) THEN   ;   jflag = 1  
    204                   ELSE                       ;   zcoef = 0._wp 
     247                  IF( zcoef > 0._wp ) THEN 
     248                     jflag = 1  
     249                  ELSE 
     250                     zcoef = 0._wp 
    205251                  ENDIF 
    206                   IF( jk1 > nn_wdit )  zcoef = 0._wp 
    207                   IF( zflxu1(ji,  jj) > 0._wp )  zwdlmtu(ji  ,jj) = zcoef 
    208                   IF( zflxu1(ji-1,jj) < 0._wp )  zwdlmtu(ji-1,jj) = zcoef 
    209                   IF( zflxv1(ji,  jj) > 0._wp )  zwdlmtv(ji  ,jj) = zcoef 
    210                   IF( zflxv1(ji,jj-1) < 0._wp )  zwdlmtv(ji,jj-1) = zcoef 
    211                ENDIF 
     252                  IF(jk1 > nn_wdit) zcoef = 0._wp 
     253                  IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef 
     254                  IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
     255                  IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
     256                  IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 
     257               END IF 
    212258            END DO ! ji loop 
    213259         END DO  ! jj loop 
    214          ! 
     260 
    215261         CALL lbc_lnk( zwdlmtu, 'U', 1. ) 
    216262         CALL lbc_lnk( zwdlmtv, 'V', 1. ) 
    217          ! 
    218          IF(lk_mpp)   CALL mpp_max(jflag)   !max over the global domain 
    219          ! 
    220          IF(jflag == 0)   EXIT 
    221          !  
     263 
     264         IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain 
     265 
     266         IF(jflag == 0) EXIT 
     267           
    222268      END DO  ! jk1 loop 
    223269        
    224270      DO jk = 1, jpkm1 
    225          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:)  
    226          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:,:)  
    227       END DO 
    228  
    229 !!gm ==> Andrew  : the lbclnk below is useless since above lbclnk is applied on zwdlmtu/v 
    230 !!                                             and un, vn always with lbclnk 
     271         un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :)  
     272         vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :)  
     273      END DO 
     274 
    231275      CALL lbc_lnk( un, 'U', -1. ) 
    232276      CALL lbc_lnk( vn, 'V', -1. ) 
    233 !!gm end 
    234       ! 
    235       un_b(:,:) = un_b(:,:) * zwdlmtu(:,:) 
    236       vn_b(:,:) = vn_b(:,:) * zwdlmtv(:,:) 
    237 !!gm ==> Andrew   : probably same as above 
     277        ! 
     278      un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 
     279      vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 
    238280      CALL lbc_lnk( un_b, 'U', -1. ) 
    239281      CALL lbc_lnk( vn_b, 'V', -1. ) 
    240 !!gm end 
    241282        
    242283      IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 
     
    245286      !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    246287      ! 
     288      ! 
     289      ! 
     290      ! 
     291      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt') 
    247292      ! 
    248293      IF( ln_timing )   CALL timing_stop('wad_lmt') 
     
    276321      !!---------------------------------------------------------------------- 
    277322      ! 
    278       IF( ln_timing )  CALL timing_start('wad_lmt_bt') 
    279       !        
    280       !IF(lwp) WRITE(numout,*) 
    281       !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting' 
    282         
     323      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt') 
     324      ! 
    283325      jflag  = 0 
    284326      zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes 
     
    296338      DO jj = 2, jpj 
    297339         DO ji = 2, jpi  
    298             ! 
    299             IF( tmask(ji,jj,1) < 0.5_wp )  CYCLE   ! we don't care about land cells 
    300             IF( ht_wd(ji,jj)   > zdepwd )   CYCLE   ! and cells which are unlikely to dry 
    301             ! 
    302             zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj) , 0._wp )  & 
    303                &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,jj-1) , 0._wp )  
    304             zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj) , 0._wp )  & 
    305                &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,jj-1) , 0._wp )  
    306  
    307             zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
     340 
     341           IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells 
     342           IF( ht_0(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry 
     343 
     344            zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + & 
     345                         & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp)  
     346            zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + & 
     347                         & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp)  
     348 
     349            zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
    308350            IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary 
    309                sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
    310                IF( zflxu(ji,  jj) > 0._wp )  zwdlmtu(ji  ,jj) = 0._wp 
    311                IF( zflxu(ji-1,jj) < 0._wp )  zwdlmtu(ji-1,jj) = 0._wp 
    312                IF( zflxv(ji,  jj) > 0._wp )  zwdlmtv(ji  ,jj) = 0._wp 
    313                IF( zflxv(ji,jj-1) < 0._wp )  zwdlmtv(ji,jj-1) = 0._wp  
    314             ENDIF 
    315          END DO 
     351              sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 
     352              IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp 
     353              IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 
     354              IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp 
     355              IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp  
     356            END IF 
     357         ENDDO 
    316358      END DO 
    317359 
     
    319361      !! start limiter iterations  
    320362      DO jk1 = 1, nn_wdit + 1 
    321          !  
     363        
     364           
    322365         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 
    323366         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    324367         jflag = 0     ! flag indicating if any further iterations are needed 
    325          ! 
     368           
    326369         DO jj = 2, jpj 
    327370            DO ji = 2, jpi  
    328                ! 
    329                IF( tmask(ji,jj, 1 ) < 0.5_wp  )  CYCLE  
    330                IF( ht_wd(ji,jj)      > zdepwd )   CYCLE 
    331                ! 
     371         
     372               IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE  
     373               IF( ht_0(ji,jj) > zdepwd )      CYCLE 
     374         
    332375               ztmp = e1e2t(ji,jj) 
    333                ! 
    334                zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj) , 0._wp )  & 
    335                   &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,jj-1) , 0._wp )  
    336                zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj) , 0._wp )  & 
    337                   &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,jj-1) , 0._wp )  
     376 
     377               zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + & 
     378                      & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp)  
     379               zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + & 
     380                      & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp)  
    338381           
    339382               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
    340                zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 
     383               zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 
    341384           
    342385               IF(zdep1 > zdep2) THEN 
    343                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    344                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 
    345                   ! flag if the limiter has been used but stop flagging if the only 
    346                   ! changes have zeroed the coefficient since further iterations will 
    347                   ! not change anything 
    348                   IF( zcoef > 0._wp ) THEN 
    349                       jflag = 1  
    350                   ELSE 
    351                       zcoef = 0._wp 
    352                   ENDIF 
    353                   IF( jk1 > nn_wdit )  zcoef = 0._wp 
    354                   IF( zflxu1(ji,  jj) > 0._wp )  zwdlmtu(ji  ,jj) = zcoef 
    355                   IF( zflxu1(ji-1,jj) < 0._wp )  zwdlmtu(ji-1,jj) = zcoef 
    356                   IF( zflxv1(ji,  jj) > 0._wp )  zwdlmtv(ji  ,jj) = zcoef 
    357                   IF( zflxv1(ji,jj-1) < 0._wp )  zwdlmtv(ji,jj-1) = zcoef 
    358                ENDIF 
     386                 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
     387                 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 
     388                 ! flag if the limiter has been used but stop flagging if the only 
     389                 ! changes have zeroed the coefficient since further iterations will 
     390                 ! not change anything 
     391                 IF( zcoef > 0._wp ) THEN 
     392                    jflag = 1  
     393                 ELSE 
     394                    zcoef = 0._wp 
     395                 ENDIF 
     396                 IF(jk1 > nn_wdit) zcoef = 0._wp 
     397                 IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef 
     398                 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
     399                 IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
     400                 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 
     401               END IF 
    359402            END DO ! ji loop 
    360403         END DO  ! jj loop 
    361          ! 
     404 
    362405         CALL lbc_lnk( zwdlmtu, 'U', 1. ) 
    363406         CALL lbc_lnk( zwdlmtv, 'V', 1. ) 
    364          ! 
    365          IF(lk_mpp)   CALL mpp_max(jflag)   !max over the global domain 
    366          ! 
    367          IF( jflag == 0 )  EXIT 
    368          !     
     407 
     408         IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain 
     409 
     410         IF(jflag == 0) EXIT 
     411           
    369412      END DO  ! jk1 loop 
    370       !  
     413       
    371414      zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :)  
    372415      zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :)  
    373       ! 
     416 
    374417      CALL lbc_lnk( zflxu, 'U', -1. ) 
    375418      CALL lbc_lnk( zflxv, 'V', -1. ) 
    376       ! 
     419        
    377420      IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 
    378421        
     
    380423      !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    381424      ! 
    382       IF( ln_timing )  CALL timing_stop('wad_lmt_bt') 
    383       ! 
     425      ! 
     426      CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 
     427      CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 
     428      ! 
     429 
     430      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt') 
    384431   END SUBROUTINE wad_lmt_bt 
    385432 
Note: See TracChangeset for help on using the changeset viewer.