Changeset 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r6152 r7646 1 2 1 MODULE wet_dry 3 2 !!============================================================================== … … 7 6 !! only effects if wetting/drying is on (ln_wd == .true.) 8 7 !!============================================================================== 9 !! History : 10 !! NEMO 3.6 ! 2014-09 ((H.Liu) Original code 8 !! History : 3.6 ! 2014-09 ((H.Liu) Original code 11 9 !! ! will add the runoff and periodic BC case later 12 10 !!---------------------------------------------------------------------- … … 33 31 !! --------------------------------------------------------------------- 34 32 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wduflt, wdvflt !: u- and v- filter36 33 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_wd !: wetting and drying t-pt depths 35 ! (can include negative depths) 37 36 38 37 LOGICAL, PUBLIC :: ln_wd !: Wetting/drying activation switch (T:on,F:off) … … 77 76 WRITE(numout,*) 78 77 WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read' 79 WRITE(numout,*) '~~~~~~~ 78 WRITE(numout,*) '~~~~~~~~' 80 79 WRITE(numout,*) ' Namelist namwad' 81 80 WRITE(numout,*) ' Logical activation ln_wd = ', ln_wd … … 84 83 WRITE(numout,*) ' land elevation threshold rn_wdld = ', rn_wdld 85 84 WRITE(numout,*) ' Max iteration for W/D limiter nn_wdit = ', nn_wdit 86 87 85 ENDIF 86 ! 88 87 IF(ln_wd) THEN 89 ALLOCATE( wd uflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj),STAT=ierr )88 ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj), STAT=ierr ) 90 89 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 91 90 ENDIF 91 ! 92 92 END SUBROUTINE wad_init 93 93 94 94 95 SUBROUTINE wad_lmt( sshb1, sshemp, z2dt ) … … 107 108 ! 108 109 INTEGER :: ji, jj, jk, jk1 ! dummy loop indices 109 INTEGER :: zflag ! local scalar110 INTEGER :: jflag ! local scalar 110 111 REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars 111 112 REAL(wp) :: zzflxp, zzflxn ! local scalars … … 116 117 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu, zflxv ! local 2D workspace 117 118 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! local 2D workspace 118 119 119 !!---------------------------------------------------------------------- 120 120 ! … … 131 131 !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting' 132 132 133 zflag = 0133 jflag = 0 134 134 zdepwd = 50._wp !maximum depth on which that W/D could possibly happen 135 135 … … 145 145 ! Horizontal Flux in u and v direction 146 146 DO jk = 1, jpkm1 147 DO jj = 1, jpj m1148 DO ji = 1, jpi m1147 DO jj = 1, jpj 148 DO ji = 1, jpi 149 149 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 150 150 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 156 156 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 157 157 158 DO jj = 2, jpjm1 159 DO ji = 2, jpim1 160 161 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells 162 IF(bathy(ji,jj) > zdepwd) CYCLE ! and cells which will unlikely go dried out 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 163 164 164 165 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & … … 167 168 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 168 169 169 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 170 IF(zdep2 < 0._wp) THEN !add more safty, but not necessary 171 !zdep2 = 0._wp 172 sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj) 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 173 178 END IF 174 179 ENDDO … … 182 187 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 183 188 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 184 185 DO jj = 2, jpjm1 186 DO ji = 2, jpim1 189 jflag = 0 ! flag indicating if any further iterations are needed 190 191 DO jj = 2, jpj 192 DO ji = 2, jpi 187 193 188 wdmask(ji,jj) = 0 189 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 190 IF(bathy(ji,jj) > zdepwd) CYCLE 194 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 195 IF( ht_wd(ji,jj) > zdepwd ) CYCLE 191 196 192 197 ztmp = e1e2t(ji,jj) … … 198 203 199 204 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 200 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) ! this one can be moved out of the loop 201 202 IF(zdep1 > zdep2) THEN 203 zflag = 1 204 wdmask(ji, jj) = 1 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 205 209 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 206 zcoef = max(zcoef, 0._wp) 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 207 219 IF(jk1 > nn_wdit) zcoef = 0._wp 208 220 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 209 221 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 210 222 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 211 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji -1,jj) = zcoef223 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 212 224 END IF 213 225 END DO ! ji loop … … 217 229 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 218 230 219 IF(lk_mpp) CALL mpp_max(zflag) !max over the global domain 220 221 IF(zflag == 0) EXIT 222 223 zflag = 0 ! flag indicating if any further iteration is needed? 231 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 232 233 IF(jflag == 0) EXIT 234 224 235 END DO ! jk1 loop 225 236 … … 231 242 CALL lbc_lnk( un, 'U', -1. ) 232 243 CALL lbc_lnk( vn, 'V', -1. ) 233 234 IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 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!!!' 235 251 236 252 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) … … 240 256 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 241 257 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 242 !243 END 244 258 ! 259 ENDIF 260 ! 245 261 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 262 ! 246 263 END SUBROUTINE wad_lmt 264 247 265 248 266 SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt ) … … 260 278 ! 261 279 INTEGER :: ji, jj, jk, jk1 ! dummy loop indices 262 INTEGER :: zflag! local scalar280 INTEGER :: jflag ! local scalar 263 281 REAL(wp) :: z2dt 264 282 REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars … … 269 287 REAL(wp), POINTER, DIMENSION(:,:) :: zflxp, zflxn ! local 2D workspace 270 288 REAL(wp), POINTER, DIMENSION(:,:) :: zflxu1, zflxv1 ! local 2D workspace 271 272 !!---------------------------------------------------------------------- 273 ! 274 289 !!---------------------------------------------------------------------- 290 ! 275 291 IF( nn_timing == 1 ) CALL timing_start('wad_lmt_bt') 276 292 … … 284 300 !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting' 285 301 286 zflag = 0302 jflag = 0 287 303 zdepwd = 50._wp !maximum depth that ocean cells can have W/D processes 288 304 … … 291 307 zflxp(:,:) = 0._wp 292 308 zflxn(:,:) = 0._wp 293 !zflxu(:,:) = 0._wp294 !zflxv(:,:) = 0._wp295 309 296 310 zwdlmtu(:,:) = 1._wp … … 299 313 ! Horizontal Flux in u and v direction 300 314 301 !zflxu(:,:) = zflxu(:,:) * e2u(:,:) 302 !zflxv(:,:) = zflxv(:,:) * e1v(:,:) 303 304 DO jj = 2, jpjm1 305 DO ji = 2, jpim1 306 307 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells 308 IF(bathy(ji,jj) > zdepwd) CYCLE ! and cells which will unlikely go dried out 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 309 320 310 321 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & … … 313 324 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 314 325 315 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 316 IF(zdep2 < 0._wp) THEN !add more safty, but not necessary 317 !zdep2 = 0._wp 318 sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj) 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 319 333 END IF 320 334 ENDDO … … 328 342 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 329 343 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 330 331 DO jj = 2, jpjm1 332 DO ji = 2, jpim1 344 jflag = 0 ! flag indicating if any further iterations are needed 345 346 DO jj = 2, jpj 347 DO ji = 2, jpi 333 348 334 wdmask(ji,jj) = 0 335 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 336 IF(bathy(ji,jj) > zdepwd) CYCLE 349 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE 350 IF( ht_wd(ji,jj) > zdepwd ) CYCLE 337 351 338 352 ztmp = e1e2t(ji,jj) … … 344 358 345 359 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 346 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 ! this one can be moved out of the loop 347 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj) 360 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 348 361 349 362 IF(zdep1 > zdep2) THEN 350 zflag = 1351 !wdmask(ji, jj) = 1352 363 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 353 zcoef = max(zcoef, 0._wp) 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 369 jflag = 1 370 ELSE 371 zcoef = 0._wp 372 ENDIF 354 373 IF(jk1 > nn_wdit) zcoef = 0._wp 355 374 IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef 356 375 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 357 376 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 358 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji -1,jj) = zcoef377 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 359 378 END IF 360 379 END DO ! ji loop … … 364 383 CALL lbc_lnk( zwdlmtv, 'V', 1. ) 365 384 366 IF(lk_mpp) CALL mpp_max(zflag) !max over the global domain 367 368 IF(zflag == 0) EXIT 369 370 zflag = 0 ! flag indicating if any further iteration is needed? 385 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 386 387 IF(jflag == 0) EXIT 388 371 389 END DO ! jk1 loop 372 390 … … 377 395 CALL lbc_lnk( zflxv, 'V', -1. ) 378 396 379 IF( zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'397 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 380 398 381 399 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) … … 385 403 CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 386 404 CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 387 !405 ! 388 406 END IF 389 407 390 408 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 391 409 END SUBROUTINE wad_lmt_bt 410 411 !!============================================================================== 392 412 END MODULE wet_dry
Note: See TracChangeset
for help on using the changeset viewer.