- Timestamp:
- 2017-12-01T18:44:09+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r7646 r8882 11 11 12 12 !!---------------------------------------------------------------------- 13 !! wad_lmt : Compute the horizontal flux limiter and the limited velocity 14 !! when wetting and drying happens 15 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and tracers 17 USE dom_oce ! ocean space and time domain 18 USE sbc_oce, ONLY : ln_rnf ! surface boundary condition: ocean 19 USE sbcrnf ! river runoff 20 USE in_out_manager ! I/O manager 21 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 22 USE lib_mpp ! MPP library 23 USE wrk_nemo ! Memory Allocation 24 USE timing ! Timing 13 !! wad_init : initialisation of wetting and drying 14 !! wad_lmt : horizontal flux limiter and limited velocity when wetting and drying happens 15 !! wad_lmt_bt : same as wad_lmt for the barotropic stepping (dynspg_ts) 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and tracers 18 USE dom_oce ! ocean space and time domain 19 USE sbc_oce , ONLY: ln_rnf ! surface boundary condition: ocean 20 USE sbcrnf ! river runoff 21 ! 22 USE in_out_manager ! I/O manager 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 USE lib_mpp ! MPP library 25 USE timing ! Timing 25 26 26 27 IMPLICIT NONE … … 31 32 !! --------------------------------------------------------------------- 32 33 33 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_wd 35 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_wd !: wetting and drying t-pt depths 36 ! ! (can include negative depths) 36 37 37 38 LOGICAL, PUBLIC :: ln_wd !: Wetting/drying activation switch (T:on,F:off) 38 39 REAL(wp), PUBLIC :: rn_wdmin1 !: minimum water depth on dried cells 39 40 REAL(wp), PUBLIC :: rn_wdmin2 !: tolerrance of minimum water depth on dried cells 40 REAL(wp), PUBLIC :: rn_wdld !: land elevation below which wetting/drying 41 !: will be considered 41 REAL(wp), PUBLIC :: rn_wdld !: land elevation below which wetting/drying will be considered 42 42 INTEGER , PUBLIC :: nn_wdit !: maximum number of iteration for W/D limiter 43 43 … … 48 48 !! * Substitutions 49 49 # include "vectopt_loop_substitute.h90" 50 !!---------------------------------------------------------------------- 50 51 CONTAINS 51 52 … … 58 59 !! ** input : - namwad namelist 59 60 !!---------------------------------------------------------------------- 61 INTEGER :: ios, ierr ! Local integer 62 !! 60 63 NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 61 INTEGER :: ios ! Local integer output status for namelist read 62 INTEGER :: ierr ! Local integer status array allocation 63 !!---------------------------------------------------------------------- 64 65 REWIND( numnam_ref ) ! Namelist namwad in reference namelist 66 ! : Parameters for Wetting/Drying 64 !!---------------------------------------------------------------------- 65 ! 66 REWIND( numnam_ref ) ! Namelist namwad in reference namelist : Parameters for Wetting/Drying 67 67 READ ( numnam_ref, namwad, IOSTAT = ios, ERR = 905) 68 68 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.) 69 REWIND( numnam_cfg ) ! Namelist namwad in configuration namelist 70 ! : Parameters for Wetting/Drying 69 REWIND( numnam_cfg ) ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying 71 70 READ ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906) 72 71 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. ) 73 72 IF(lwm) WRITE ( numond, namwad ) 74 73 ! 75 74 IF(lwp) THEN ! control print 76 75 WRITE(numout,*) … … 103 102 !! ** Action : - calculate flux limiter and W/D flag 104 103 !!---------------------------------------------------------------------- 105 REAL(wp), DIMENSION(:,:), INTENT(inout) :: sshb1 106 REAL(wp), DIMENSION(:,:), INTENT(in ):: sshemp107 REAL(wp) , INTENT(in) ::z2dt104 REAL(wp), DIMENSION(:,:), INTENT(inout) :: sshb1 !!gm DOCTOR names: should start with p ! 105 REAL(wp), DIMENSION(:,:), INTENT(in ) :: sshemp 106 REAL(wp) , INTENT(in ) :: z2dt 108 107 ! 109 108 INTEGER :: ji, jj, jk, jk1 ! dummy loop indices … … 113 112 REAL(wp) :: zdepwd ! local scalar, always wet cell depth 114 113 REAL(wp) :: ztmp ! local scalars 115 REAL(wp), POINTER, DIMENSION(:,:) :: zwdlmtu, zwdlmtv !: W/D flux limiters 116 REAL(wp), POINTER, DIMENSION(:,:) :: zflxp, zflxn ! local 2D workspace 117 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu, zflxv ! local 2D workspace 118 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! local 2D workspace 119 !!---------------------------------------------------------------------- 120 ! 121 122 IF( nn_timing == 1 ) CALL timing_start('wad_lmt') 123 124 IF(ln_wd) THEN 125 126 CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 127 CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv) 128 ! 129 130 !IF(lwp) WRITE(numout,*) 131 !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting' 132 133 jflag = 0 134 zdepwd = 50._wp !maximum depth on which that W/D could possibly happen 135 136 137 zflxp(:,:) = 0._wp 138 zflxn(:,:) = 0._wp 139 zflxu(:,:) = 0._wp 140 zflxv(:,:) = 0._wp 141 142 zwdlmtu(:,:) = 1._wp 143 zwdlmtv(:,:) = 1._wp 144 145 ! Horizontal Flux in u and v direction 146 DO jk = 1, jpkm1 147 DO jj = 1, jpj 148 DO ji = 1, jpi 149 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 150 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 151 END DO 152 END DO 153 END DO 154 155 zflxu(:,:) = zflxu(:,:) * e2u(:,:) 156 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 157 158 wdmask(:,:) = 1 159 DO jj = 2, jpj 160 DO ji = 2, jpi 161 162 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE ! we don't care about land cells 163 IF( ht_wd(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry 164 165 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & 166 & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp) 167 zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) + & 168 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 169 170 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 171 IF(zdep2 .le. 0._wp) THEN !add more safty, but not necessary 172 sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 173 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 174 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 175 IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp 176 IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 177 wdmask(ji,jj) = 0._wp 178 END IF 179 ENDDO 180 END DO 181 182 183 !! start limiter iterations 184 DO jk1 = 1, nn_wdit + 1 185 186 187 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 188 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 189 jflag = 0 ! flag indicating if any further iterations are needed 190 191 DO jj = 2, jpj 192 DO ji = 2, jpi 193 194 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 195 IF( ht_wd(ji,jj) > zdepwd ) CYCLE 196 197 ztmp = e1e2t(ji,jj) 198 199 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + & 200 & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp) 201 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 202 & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp) 203 204 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 205 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 206 207 IF( zdep1 > zdep2 ) THEN 208 wdmask(ji, jj) = 0 209 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 210 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 211 ! flag if the limiter has been used but stop flagging if the only 212 ! changes have zeroed the coefficient since further iterations will 213 ! not change anything 214 IF( zcoef > 0._wp ) THEN 215 jflag = 1 216 ELSE 217 zcoef = 0._wp 218 ENDIF 219 IF(jk1 > nn_wdit) zcoef = 0._wp 220 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 221 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 222 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 223 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 224 END IF 225 END DO ! ji loop 226 END DO ! jj loop 227 228 CALL lbc_lnk( zwdlmtu, 'U', 1. ) 229 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 230 231 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 232 233 IF(jflag == 0) EXIT 234 235 END DO ! jk1 loop 236 237 DO jk = 1, jpkm1 238 un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 239 vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 240 END DO 241 242 CALL lbc_lnk( un, 'U', -1. ) 243 CALL lbc_lnk( vn, 'V', -1. ) 244 ! 245 un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 246 vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 247 CALL lbc_lnk( un_b, 'U', -1. ) 248 CALL lbc_lnk( vn_b, 'V', -1. ) 249 250 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 251 252 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 253 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 254 ! 255 ! 256 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 257 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 258 ! 259 ENDIF 260 ! 261 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 114 REAL(wp), DIMENSION(jpi,jpj) :: zwdlmtu, zwdlmtv ! W/D flux limiters 115 REAL(wp), DIMENSION(jpi,jpj) :: zflxp , zflxn ! local 2D workspace 116 REAL(wp), DIMENSION(jpi,jpj) :: zflxu , zflxv ! local 2D workspace 117 REAL(wp), DIMENSION(jpi,jpj) :: zflxu1 , zflxv1 ! local 2D workspace 118 !!---------------------------------------------------------------------- 119 ! 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 ! 125 jflag = 0 126 zdepwd = 50._wp !maximum depth on which that W/D could possibly happen 127 ! 128 zflxp(:,:) = 0._wp 129 zflxn(:,:) = 0._wp 130 zflxu(:,:) = 0._wp 131 zflxv(:,:) = 0._wp 132 ! 133 zwdlmtu(:,:) = 1._wp 134 zwdlmtv(:,:) = 1._wp 135 ! 136 ! Horizontal Flux in u and v direction 137 DO jk = 1, jpkm1 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 141 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 142 END DO 143 END DO 144 END DO 145 ! 146 zflxu(:,:) = zflxu(:,:) * e2u(:,:) 147 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 148 ! 149 wdmask(:,:) = 1 150 DO jj = 2, jpj 151 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 168 wdmask(ji,jj) = 0._wp 169 ENDIF 170 END DO 171 END DO 172 !! 173 !! start limiter iterations 174 DO jk1 = 1, nn_wdit + 1 175 ! 176 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 177 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 178 jflag = 0 ! flag indicating if any further iterations are needed 179 ! 180 DO jj = 2, jpj 181 DO ji = 2, jpi 182 ! 183 IF( tmask(ji,jj,1) < 0.5_wp ) CYCLE 184 IF( ht_wd(ji,jj) > zdepwd ) CYCLE 185 ! 186 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 ! 193 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 194 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 195 ! 196 IF( zdep1 > zdep2 ) THEN 197 wdmask(ji, jj) = 0 198 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 199 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 200 ! flag if the limiter has been used but stop flagging if the only 201 ! changes have zeroed the coefficient since further iterations will 202 ! not change anything 203 IF( zcoef > 0._wp ) THEN ; jflag = 1 204 ELSE ; zcoef = 0._wp 205 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 212 END DO ! ji loop 213 END DO ! jj loop 214 ! 215 CALL lbc_lnk( zwdlmtu, 'U', 1. ) 216 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 ! 222 END DO ! jk1 loop 223 224 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 231 CALL lbc_lnk( un, 'U', -1. ) 232 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 238 CALL lbc_lnk( un_b, 'U', -1. ) 239 CALL lbc_lnk( vn_b, 'V', -1. ) 240 !!gm end 241 242 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 243 244 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 245 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 246 ! 247 ! 248 IF( ln_timing ) CALL timing_stop('wad_lmt') 262 249 ! 263 250 END SUBROUTINE wad_lmt … … 284 271 REAL(wp) :: zdepwd ! local scalar, always wet cell depth 285 272 REAL(wp) :: ztmp ! local scalars 286 REAL(wp), POINTER, DIMENSION(:,:) :: zwdlmtu, zwdlmtv !: W/D flux limiters 287 REAL(wp), POINTER, DIMENSION(:,:) :: zflxp, zflxn ! local 2D workspace 288 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! local 2D workspace 289 !!---------------------------------------------------------------------- 290 ! 291 IF( nn_timing == 1 ) CALL timing_start('wad_lmt_bt') 292 293 IF(ln_wd) THEN 294 295 CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 296 CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv) 297 ! 298 299 !IF(lwp) WRITE(numout,*) 300 !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting' 301 302 jflag = 0 303 zdepwd = 50._wp !maximum depth that ocean cells can have W/D processes 304 305 z2dt = rdtbt 306 307 zflxp(:,:) = 0._wp 308 zflxn(:,:) = 0._wp 309 310 zwdlmtu(:,:) = 1._wp 311 zwdlmtv(:,:) = 1._wp 312 313 ! Horizontal Flux in u and v direction 314 315 DO jj = 2, jpj 316 DO ji = 2, jpi 317 318 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells 319 IF( ht_wd(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry 320 321 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & 322 & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp) 323 zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) + & 324 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 325 326 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 327 IF(zdep2 .le. 0._wp) THEN !add more safety, but not necessary 328 sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 329 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 330 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 331 IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp 332 IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 333 END IF 334 ENDDO 335 END DO 273 REAL(wp), DIMENSION(jpi,jpj) :: zwdlmtu, zwdlmtv !: W/D flux limiters 274 REAL(wp), DIMENSION(jpi,jpj) :: zflxp, zflxn ! local 2D workspace 275 REAL(wp), DIMENSION(jpi,jpj) :: zflxu1, zflxv1 ! local 2D workspace 276 !!---------------------------------------------------------------------- 277 ! 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 283 jflag = 0 284 zdepwd = 50._wp !maximum depth that ocean cells can have W/D processes 285 286 z2dt = rdtbt 287 288 zflxp(:,:) = 0._wp 289 zflxn(:,:) = 0._wp 290 291 zwdlmtu(:,:) = 1._wp 292 zwdlmtv(:,:) = 1._wp 293 294 ! Horizontal Flux in u and v direction 295 296 DO jj = 2, jpj 297 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 308 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 316 END DO 336 317 337 318 338 !! start limiter iterations 339 DO jk1 = 1, nn_wdit + 1 340 319 !! start limiter iterations 320 DO jk1 = 1, nn_wdit + 1 321 ! 322 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 323 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 324 jflag = 0 ! flag indicating if any further iterations are needed 325 ! 326 DO jj = 2, jpj 327 DO ji = 2, jpi 328 ! 329 IF( tmask(ji,jj, 1 ) < 0.5_wp ) CYCLE 330 IF( ht_wd(ji,jj) > zdepwd ) CYCLE 331 ! 332 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 ) 341 338 342 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 343 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 344 jflag = 0 ! flag indicating if any further iterations are needed 339 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 340 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 345 341 346 DO jj = 2, jpj 347 DO ji = 2, jpi 348 349 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE 350 IF( ht_wd(ji,jj) > zdepwd ) CYCLE 351 352 ztmp = e1e2t(ji,jj) 353 354 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + & 355 & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp) 356 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 357 & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp) 358 359 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 360 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 361 362 IF(zdep1 > zdep2) THEN 363 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 364 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 365 ! flag if the limiter has been used but stop flagging if the only 366 ! changes have zeroed the coefficient since further iterations will 367 ! not change anything 368 IF( zcoef > 0._wp ) THEN 342 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 369 349 jflag = 1 370 350 ELSE 371 351 zcoef = 0._wp 372 ENDIF 373 IF(jk1 > nn_wdit) zcoef = 0._wp 374 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 375 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 376 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 377 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 378 END IF 379 END DO ! ji loop 380 END DO ! jj loop 381 382 CALL lbc_lnk( zwdlmtu, 'U', 1. ) 383 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 384 385 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 386 387 IF(jflag == 0) EXIT 388 389 END DO ! jk1 loop 390 391 zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 392 zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 393 394 CALL lbc_lnk( zflxu, 'U', -1. ) 395 CALL lbc_lnk( zflxv, 'V', -1. ) 396 397 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 398 399 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 400 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 401 ! 402 ! 403 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 404 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 405 ! 406 END IF 407 408 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 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 359 END DO ! ji loop 360 END DO ! jj loop 361 ! 362 CALL lbc_lnk( zwdlmtu, 'U', 1. ) 363 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 ! 369 END DO ! jk1 loop 370 ! 371 zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 372 zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 373 ! 374 CALL lbc_lnk( zflxu, 'U', -1. ) 375 CALL lbc_lnk( zflxv, 'V', -1. ) 376 ! 377 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 378 379 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 380 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 381 ! 382 IF( ln_timing ) CALL timing_stop('wad_lmt') 383 ! 409 384 END SUBROUTINE wad_lmt_bt 410 385
Note: See TracChangeset
for help on using the changeset viewer.
