- Timestamp:
- 2017-12-19T09:26:25+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r9092 r9124 26 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 27 USE lib_mpp ! MPP library 28 USE timing ! Timing28 USE timing ! timing of the main modules 29 29 30 30 IMPLICIT NONE … … 71 71 !! ** input : - namwad namelist 72 72 !!---------------------------------------------------------------------- 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 73 INTEGER :: ios, ierr ! Local integer 74 !! 75 NAMELIST/namwad/ ln_wd_il, ln_wd_dl , rn_wdmin0, rn_wdmin1, rn_wdmin2, rn_wdld, & 76 & nn_wdit , ln_wd_dl_bc, ln_wd_dl_rmp 78 77 !!---------------------------------------------------------------------- 79 78 ! … … 100 99 WRITE(numout,*) ' T => baroclinic u,v=0 at dry pts: ln_wd_dl_bc = ', ln_wd_dl_bc 101 100 WRITE(numout,*) ' use a ramp for rwd limiter: ln_wd_dl_rwd_rmp = ', ln_wd_dl_rmp 102 103 101 ENDIF 104 102 IF( .NOT. ln_read_cfg ) THEN 105 103 IF(lwp) WRITE(numout,*) ' No configuration file so seting ssh_ref to zero ' 106 ssh_ref=0. 0104 ssh_ref=0._wp 107 105 ENDIF 108 106 109 r_rn_wdmin1 =1/rn_wdmin1107 r_rn_wdmin1 = 1 / rn_wdmin1 110 108 ll_wd = .FALSE. 111 IF( ln_wd_il .OR. ln_wd_dl) THEN109 IF( ln_wd_il .OR. ln_wd_dl ) THEN 112 110 ll_wd = .TRUE. 113 111 ALLOCATE( wdmask(jpi,jpj), STAT=ierr ) … … 144 142 REAL(wp), DIMENSION(jpi,jpj) :: zflxu1 , zflxv1 ! local 2D workspace 145 143 !!---------------------------------------------------------------------- 146 ! 147 IF( nn_timing == 1 ) CALL timing_start('wad_lmt') 148 ! 149 144 IF( ln_timing ) CALL timing_start('wad_lmt') ! 145 ! 150 146 DO jk = 1, jpkm1 151 147 un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:) … … 153 149 END DO 154 150 jflag = 0 155 zdepwd = 50._wp !maximum depth on which that W/D could possibly happen156 151 zdepwd = 50._wp ! maximum depth on which that W/D could possibly happen 152 ! 157 153 zflxp(:,:) = 0._wp 158 154 zflxn(:,:) = 0._wp 159 155 zflxu(:,:) = 0._wp 160 156 zflxv(:,:) = 0._wp 161 162 zwdlmtu(:,:) = 1._wp 163 zwdlmtv(:,:) = 1._wp 164 165 ! Horizontal Flux in u and v direction 166 DO jk = 1, jpkm1 167 DO jj = 1, jpj 168 DO ji = 1, jpi 169 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 170 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 171 END DO 172 END DO 173 END DO 174 157 ! 158 zwdlmtu(:,:) = 1._wp 159 zwdlmtv(:,:) = 1._wp 160 ! 161 DO jk = 1, jpkm1 ! Horizontal Flux in u and v direction 162 zflxu(:,:) = zflxu(:,:) + e3u_n(:,:,jk) * un(:,:,jk) * umask(:,:,jk) 163 zflxv(:,:) = zflxv(:,:) + e3v_n(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk) 164 END DO 175 165 zflxu(:,:) = zflxu(:,:) * e2u(:,:) 176 166 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 177 178 wdmask(:,:) = 1 167 ! 168 wdmask(:,:) = 1._wp 179 169 DO jj = 2, jpj 180 170 DO ji = 2, jpi 181 182 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE! we don't care about land cells183 IF( ht_0(ji,jj) - ssh_ref > zdepwd ) CYCLE! and cells which are unlikely to dry184 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 171 ! 172 IF( tmask(ji,jj,1) < 0.5_wp ) CYCLE ! we don't care about land cells 173 IF( ht_0(ji,jj) - ssh_ref > zdepwd ) CYCLE ! and cells which are unlikely to dry 174 ! 175 zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj ) , 0._wp ) & 176 & + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji, jj-1) , 0._wp ) 177 zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj ) , 0._wp ) & 178 & + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji, jj-1) , 0._wp ) 179 ! 190 180 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 191 IF( zdep2 .le. 0._wp) THEN !add more safty, but not necessary181 IF( zdep2 <= 0._wp ) THEN ! add more safty, but not necessary 192 182 sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 193 183 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp … … 197 187 wdmask(ji,jj) = 0._wp 198 188 END IF 199 ENDDO 200 END DO 201 202 203 ! HPG limiter from jholt 189 END DO 190 END DO 191 ! 192 ! ! HPG limiter from jholt 204 193 wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp) 205 !jth assume don't need a lbc_lnk here194 !jth assume don't need a lbc_lnk here 206 195 DO jj = 1, jpjm1 207 196 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))197 wdrampu(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji+1,jj) ) 198 wdrampv(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji,jj+1) ) 210 199 END DO 211 200 END DO 212 ! end HPG limiter 213 214 215 216 !! start limiter iterations 217 DO jk1 = 1, nn_wdit + 1 218 219 201 ! ! end HPG limiter 202 ! 203 ! 204 DO jk1 = 1, nn_wdit + 1 !== start limiter iterations ==! 205 ! 220 206 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 221 207 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 222 208 jflag = 0 ! flag indicating if any further iterations are needed 223 209 ! 224 210 DO jj = 2, jpj 225 211 DO ji = 2, jpi 226 227 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 228 IF( ht_0(ji,jj) > zdepwd ) CYCLE 229 212 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 213 IF( ht_0(ji,jj) > zdepwd ) CYCLE 214 ! 230 215 ztmp = e1e2t(ji,jj) 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 216 ! 217 zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj ) , 0._wp) & 218 & + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji, jj-1) , 0._wp) 219 zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj ) , 0._wp) & 220 & + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji, jj-1) , 0._wp) 221 ! 237 222 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 238 223 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 239 224 ! 240 225 IF( zdep1 > zdep2 ) THEN 241 wdmask(ji, jj) = 0 226 wdmask(ji, jj) = 0._wp 242 227 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 243 228 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) … … 245 230 ! changes have zeroed the coefficient since further iterations will 246 231 ! not change anything 247 IF( zcoef > 0._wp ) THEN 248 jflag = 1 249 ELSE 250 zcoef = 0._wp 232 IF( zcoef > 0._wp ) THEN ; jflag = 1 233 ELSE ; zcoef = 0._wp 251 234 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 258 END DO ! ji loop 259 END DO ! jj loop 260 235 IF( jk1 > nn_wdit ) zcoef = 0._wp 236 IF( zflxu1(ji ,jj ) > 0._wp ) zwdlmtu(ji ,jj ) = zcoef 237 IF( zflxu1(ji-1,jj ) < 0._wp ) zwdlmtu(ji-1,jj ) = zcoef 238 IF( zflxv1(ji ,jj ) > 0._wp ) zwdlmtv(ji ,jj ) = zcoef 239 IF( zflxv1(ji ,jj-1) < 0._wp ) zwdlmtv(ji ,jj-1) = zcoef 240 ENDIF 241 END DO 242 END DO 261 243 CALL lbc_lnk_multi( zwdlmtu, 'U', 1., zwdlmtv, 'V', 1. ) 262 263 IF( lk_mpp)CALL mpp_max(jflag) !max over the global domain264 265 IF( jflag == 0)EXIT266 244 ! 245 IF( lk_mpp ) CALL mpp_max(jflag) !max over the global domain 246 ! 247 IF( jflag == 0 ) EXIT 248 ! 267 249 END DO ! jk1 loop 268 250 ! 269 251 DO jk = 1, jpkm1 270 un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 271 vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 272 END DO 273 274 CALL lbc_lnk_multi( un, 'U', -1., vn, 'V', -1. ) 275 ! 252 un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:) 253 vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:,:) 254 END DO 276 255 un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 277 256 vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 257 ! 258 !!gm TO BE SUPPRESSED ? these lbc_lnk are useless since zwdlmtu and zwdlmtv are defined everywhere ! 259 CALL lbc_lnk_multi( un , 'U', -1., vn , 'V', -1. ) 278 260 CALL lbc_lnk_multi( un_b, 'U', -1., vn_b, 'V', -1. ) 279 280 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 281 261 !!gm 262 ! 263 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 264 ! 282 265 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 283 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 284 ! 285 ! 286 ! 287 ! 288 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 289 ! 290 IF( ln_timing ) CALL timing_stop('wad_lmt') 266 ! 267 IF( ln_timing ) CALL timing_stop('wad_lmt') ! 291 268 ! 292 269 END SUBROUTINE wad_lmt … … 303 280 !! ** Action : - calculate flux limiter and W/D flag 304 281 !!---------------------------------------------------------------------- 305 REAL(wp) , INTENT(in):: rdtbt ! ocean time-step index282 REAL(wp) , INTENT(in ) :: rdtbt ! ocean time-step index 306 283 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zflxu, zflxv, sshn_e, zssh_frc 307 284 ! 308 285 INTEGER :: ji, jj, jk, jk1 ! dummy loop indices 309 INTEGER :: jflag ! local scalar286 INTEGER :: jflag ! local integer 310 287 REAL(wp) :: z2dt 311 288 REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars … … 317 294 REAL(wp), DIMENSION(jpi,jpj) :: zflxu1, zflxv1 ! local 2D workspace 318 295 !!---------------------------------------------------------------------- 319 ! 320 IF( nn_timing == 1 ) CALL timing_start('wad_lmt_bt') 296 IF( ln_timing ) CALL timing_start('wad_lmt_bt') ! 321 297 ! 322 298 jflag = 0 323 zdepwd = 50._wp ! maximum depth that ocean cells can have W/D processes324 299 zdepwd = 50._wp ! maximum depth that ocean cells can have W/D processes 300 ! 325 301 z2dt = rdtbt 326 302 ! 327 303 zflxp(:,:) = 0._wp 328 304 zflxn(:,:) = 0._wp 329 330 305 zwdlmtu(:,:) = 1._wp 331 306 zwdlmtv(:,:) = 1._wp 332 333 ! Horizontal Flux in u and v direction 334 335 DO jj = 2, jpj 307 ! 308 DO jj = 2, jpj ! Horizontal Flux in u and v direction 336 309 DO ji = 2, jpi 337 338 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells339 IF( ht_0(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry340 341 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) +&342 & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp)343 zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) +&344 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp)345 310 ! 311 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells 312 IF( ht_0(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry 313 ! 314 zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj ) , 0._wp ) & 315 & + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji, jj-1) , 0._wp ) 316 zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj ) , 0._wp ) & 317 & + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji, jj-1) , 0._wp ) 318 ! 346 319 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 347 IF( zdep2 .le. 0._wp) THEN !add more safety, but not necessary320 IF( zdep2 <= 0._wp ) THEN !add more safety, but not necessary 348 321 sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 349 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 350 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 351 IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp 352 IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 353 END IF 354 ENDDO 355 END DO 356 357 358 !! start limiter iterations 359 DO jk1 = 1, nn_wdit + 1 360 361 322 IF( zflxu(ji ,jj ) > 0._wp) zwdlmtu(ji ,jj ) = 0._wp 323 IF( zflxu(ji-1,jj ) < 0._wp) zwdlmtu(ji-1,jj ) = 0._wp 324 IF( zflxv(ji ,jj ) > 0._wp) zwdlmtv(ji ,jj ) = 0._wp 325 IF( zflxv(ji ,jj-1) < 0._wp) zwdlmtv(ji ,jj-1) = 0._wp 326 ENDIF 327 END DO 328 END DO 329 ! 330 DO jk1 = 1, nn_wdit + 1 !! start limiter iterations 331 ! 362 332 zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 363 333 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 364 334 jflag = 0 ! flag indicating if any further iterations are needed 365 335 ! 366 336 DO jj = 2, jpj 367 337 DO ji = 2, jpi 368 369 IF( tmask(ji, jj, 1 ) < 0.5_wp )CYCLE370 IF( ht_0(ji,jj) > zdepwd )CYCLE371 338 ! 339 IF( tmask(ji, jj, 1 ) < 0.5_wp ) CYCLE 340 IF( ht_0(ji,jj) > zdepwd ) CYCLE 341 ! 372 342 ztmp = e1e2t(ji,jj) 373 374 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) +&375 &max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp)376 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) +&377 &min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp)343 ! 344 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) & 345 & + max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp) 346 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) & 347 & + min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp) 378 348 379 349 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp … … 399 369 END DO ! ji loop 400 370 END DO ! jj loop 401 371 ! 402 372 CALL lbc_lnk_multi( zwdlmtu, 'U', 1., zwdlmtv, 'V', 1. ) 403 373 ! 404 374 IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain 405 406 IF(jflag == 0) EXIT407 375 ! 376 IF(jflag == 0) EXIT 377 ! 408 378 END DO ! jk1 loop 409 379 ! 410 380 zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 411 381 zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 412 382 ! 383 !!gm THIS lbc_lnk is useless since it is already done at the end of the jk1-loop 413 384 CALL lbc_lnk_multi( zflxu, 'U', -1., zflxv, 'V', -1. ) 414 415 IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 416 385 !!gm end 386 ! 387 IF( jflag == 1 .AND. lwp ) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 388 ! 417 389 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 418 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 419 ! 420 ! 421 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 390 ! 391 IF( ln_timing ) CALL timing_stop('wad_lmt_bt') ! 392 ! 422 393 END SUBROUTINE wad_lmt_bt 423 394
Note: See TracChangeset
for help on using the changeset viewer.
