Changeset 8865 for branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN
- Timestamp:
- 2017-12-01T05:41:32+01:00 (7 years ago)
- Location:
- branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 4 edited
-
dynhpg.F90 (modified) (14 diffs)
-
dynnxt.F90 (modified) (2 diffs)
-
dynspg_ts.F90 (modified) (26 diffs)
-
wet_dry.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r7761 r8865 455 455 DO ji = 2, jpim1 456 456 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &458 & MAX( sshn(ji,jj) + ht_ wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &457 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 458 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 459 459 & > rn_wdmin1 + rn_wdmin2 460 460 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 461 461 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 462 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )462 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 463 463 464 464 IF(ll_tmp1) THEN … … 466 466 ELSE IF(ll_tmp2) THEN 467 467 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 468 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_ wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &468 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 469 469 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 470 470 ELSE … … 473 473 474 474 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 475 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &476 & MAX( sshn(ji,jj) + ht_ wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &475 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 476 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 477 477 & > rn_wdmin1 + rn_wdmin2 478 478 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 479 479 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 480 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )480 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 481 481 482 482 IF(ll_tmp1) THEN … … 484 484 ELSE IF(ll_tmp2) THEN 485 485 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 486 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_ wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &486 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 487 487 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 488 488 ELSE … … 707 707 DO ji = 2, jpim1 708 708 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 709 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &710 & MAX( sshn(ji,jj) + ht_ wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &709 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 710 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 711 711 & > rn_wdmin1 + rn_wdmin2 712 712 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 713 713 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 714 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )714 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 715 715 716 716 IF(ll_tmp1) THEN … … 718 718 ELSE IF(ll_tmp2) THEN 719 719 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 720 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_ wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &720 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 721 721 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 722 722 ELSE … … 725 725 726 726 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 727 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &728 & MAX( sshn(ji,jj) + ht_ wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &727 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 728 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 729 729 & > rn_wdmin1 + rn_wdmin2 730 730 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 731 731 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 732 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )732 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 733 733 734 734 IF(ll_tmp1) THEN … … 736 736 ELSE IF(ll_tmp2) THEN 737 737 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 738 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_ wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &738 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 739 739 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 740 740 ELSE … … 1006 1006 DO ji = 2, jpim1 1007 1007 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1008 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &1009 & MAX( sshn(ji,jj) + ht_ wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &1008 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1009 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1010 1010 & > rn_wdmin1 + rn_wdmin2 1011 1011 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( & 1012 1012 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1013 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )1013 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1014 1014 1015 1015 IF(ll_tmp1) THEN … … 1017 1017 ELSE IF(ll_tmp2) THEN 1018 1018 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1019 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_ wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &1019 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 1020 1020 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1021 1022 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1021 1023 ELSE 1022 1024 zcpx(ji,jj) = 0._wp … … 1024 1026 1025 1027 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1026 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &1027 & MAX( sshn(ji,jj) + ht_ wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &1028 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1029 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1028 1030 & > rn_wdmin1 + rn_wdmin2 1029 1031 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( & 1030 1032 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1031 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )1033 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1032 1034 1033 1035 IF(ll_tmp1) THEN … … 1035 1037 ELSE IF(ll_tmp2) THEN 1036 1038 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1037 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_ wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &1039 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1038 1040 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1041 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1042 1039 1043 ELSE 1040 1044 zcpy(ji,jj) = 0._wp … … 1228 1232 ENDIF 1229 1233 IF( ln_wd ) THEN 1230 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1231 zdpdx2 = zdpdx2 * zcpx(ji,jj) 1234 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 1235 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1232 1236 ENDIF 1233 1237 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) … … 1287 1291 ENDIF 1288 1292 IF( ln_wd ) THEN 1289 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1290 zdpdy2 = zdpdy2 * zcpy(ji,jj) 1293 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1294 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1291 1295 ENDIF 1292 1296 -
branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r7753 r8865 28 28 USE dom_oce ! ocean space and time domain 29 29 USE sbc_oce ! Surface boundary condition: ocean fields 30 USE sbcrnf ! river runoffs 30 31 USE phycst ! physical constants 31 32 USE dynadv ! dynamics: vector invariant versus flux form … … 220 221 IF ( .NOT. ln_isf ) THEN ! if no ice shelf melting 221 222 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) & 222 & - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 223 & ) * tmask(:,:,1) 224 IF( ln_rnf_depth ) THEN 225 DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too, not sure for ice shelf case yet 226 DO jj = 1, jpj 227 DO ji = 1, jpi 228 IF( mikt(ji,jj) <= jk .and. jk <= nk_rnf(ji,jj) ) THEN 229 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * (rnf_b(ji,jj) - rnf(ji,jj))*(e3t_n(ji,jj,jk)/h_rnf(ji,jj) )*tmask(ji,jj,jk) 230 ENDIF 231 ENDDO 232 ENDDO 233 ENDDO 234 ELSE 235 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * (rnf_b(ji,jj) - rnf(ji,jj))*tmask(ji,jj,1) 236 ENDIF 223 237 ELSE ! if ice shelf melting 224 238 DO jj = 1, jpj -
branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7831 r8865 1 1 MODULE dynspg_ts 2 3 !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out ! 4 2 5 !!====================================================================== 3 6 !! *** MODULE dynspg_ts *** … … 150 153 REAL(wp) :: zhura, zhvra ! - - 151 154 REAL(wp) :: za0, za1, za2, za3 ! - - 155 REAL(wp) :: zwdramp ! local scalar - only used if ln_rwd = .True. 156 157 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point 158 152 159 ! 153 160 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e … … 158 165 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 159 166 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 167 REAL(wp), POINTER, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 168 REAL(wp), POINTER, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask 160 169 !!---------------------------------------------------------------------- 161 170 ! … … 170 179 CALL wrk_alloc( jpi,jpj, zhf ) 171 180 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 181 IF( ln_rwd ) CALL wrk_alloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2) 182 183 IF ( ln_wd_diag ) THEN 184 iwdg = jn_wd_i ; jwdg = jn_wd_j ; kwdg = jn_wd_k 185 WRITE(numout,*) 'kt, iwdg, jwdg, kwdg = ', kt, iwdg, jwdg, kwdg 186 END IF 187 172 188 ! 173 189 zmdi=1.e+20 ! missing data indicator for masking … … 178 194 z1_2 = 0.5_wp 179 195 zraur = 1._wp / rau0 196 zwdramp = 1._wp / rn_wdmin1 ! simplest ramp 197 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 180 198 ! ! reciprocal of baroclinic time step 181 199 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt … … 403 421 DO ji = 2, jpim1 404 422 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 405 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &406 & MAX( sshn(ji,jj) + ht_ wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &423 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 424 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 407 425 & > rn_wdmin1 + rn_wdmin2 408 426 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 409 427 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 410 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )428 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 411 429 412 430 IF(ll_tmp1) THEN … … 414 432 ELSE IF(ll_tmp2) THEN 415 433 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 416 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_ wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &434 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 417 435 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 436 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 437 418 438 ELSE 419 439 zcpx(ji,jj) = 0._wp … … 421 441 422 442 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 423 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &424 & MAX( sshn(ji,jj) + ht_ wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &443 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 444 & MAX( sshn(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 425 445 & > rn_wdmin1 + rn_wdmin2 426 446 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 427 447 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 428 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )448 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 429 449 430 450 IF(ll_tmp1) THEN … … 432 452 ELSE IF(ll_tmp2) THEN 433 453 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 434 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_ wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &454 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 435 455 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 456 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 457 436 458 ELSE 437 459 zcpy(ji,jj) = 0._wp … … 443 465 DO ji = 2, jpim1 444 466 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 445 & * r1_e1u(ji,jj) * zcpx(ji,jj) 467 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 446 468 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 447 & * r1_e2v(ji,jj) * zcpy(ji,jj) 469 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 470 448 471 END DO 449 472 END DO … … 491 514 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 492 515 IF( ln_wd ) THEN 493 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 494 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 516 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) * wdrampu(ji,jj) 517 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) * wdrampv(ji,jj) 495 518 ELSE 496 519 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) … … 627 650 vn_adv(:,:) = 0._wp 628 651 ! ! ==================== ! 652 653 IF (ln_rwd) THEN 654 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 655 zvwdmask(:,:) = 0._wp ! 656 zuwdav2(:,:) = 0._wp 657 zvwdav2(:,:) = 0._wp 658 END IF 659 660 629 661 DO jn = 1, icycle ! sub-time-step loop ! 630 662 ! ! ==================== ! … … 654 686 ! Extrapolate Sea Level at step jit+0.5: 655 687 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 656 ! 688 689 ! set wetting & drying mask at tracer points for this barotropic sub-step 690 IF ( ln_rwd ) THEN 691 692 IF ( ln_rwd_rmp ) THEN 693 DO jj = 1, jpj 694 DO ji = 1, jpi ! vector opt. 695 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 696 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 697 ztwdmask(ji,jj) = 1._wp 698 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 699 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )/rn_wdmin1)) ) 700 ELSE 701 ztwdmask(ji,jj) = 0._wp 702 END IF 703 END DO 704 END DO 705 ELSE 706 DO jj = 1, jpj 707 DO ji = 1, jpi ! vector opt. 708 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 709 ztwdmask(ji,jj) = 1._wp 710 ELSE 711 ztwdmask(ji,jj) = 0._wp 712 END IF 713 END DO 714 END DO 715 END IF 716 717 IF ( ln_wd_diag ) WRITE(numout,*) 'kt, jn = ', kt, jn 718 IF ( ln_wd_diag ) WRITE(numout, *) 'zsshp2_e: (i,j), (i+1,j), (i,j+1) = ', zsshp2_e(iwdg,jwdg), zsshp2_e(iwdg+1,jwdg), zsshp2_e(iwdg,jwdg+1) 719 IF ( ln_wd_diag ) WRITE(numout, *) 'ht_0: (i,j), (i+1,j), (i,j+1) = ', ht_0(iwdg,jwdg), ht_0(iwdg+1,jwdg), (iwdg,jwdg+1) 720 IF ( ln_wd_diag ) WRITE(numout, *) 'ztwdmask: (i,j), (i+1,j), (i,j+1) = ', ztwdmask(iwdg,jwdg), ztwdmask(iwdg+1,jwdg), ztwdmask(iwdg,jwdg+1) 721 END IF 722 723 657 724 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 658 725 DO ji = 2, fs_jpim1 ! Vector opt. … … 707 774 #endif 708 775 IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 709 ! 776 777 IF ( ln_rwd ) THEN 778 779 IF ( ln_wd_diag ) THEN 780 WRITE(numout, *) 'zwx: (i,j), (i+1,j) = ', zwx(iwdg,jwdg), zwx(iwdg+1,jwdg) 781 WRITE(numout, *) 'zwy: (i,j), (i,j+1) = ', zwy(iwdg,jwdg), zwx(iwdg,jwdg+1) 782 END IF 783 784 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 785 786 DO jj = 1, jpjm1 787 DO ji = 1, jpim1 788 IF ( zwx(ji,jj) > 0.0 ) THEN 789 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 790 ELSE 791 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 792 END IF 793 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 794 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 795 796 IF ( zwy(ji,jj) > 0.0 ) THEN 797 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 798 ELSE 799 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 800 END IF 801 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 802 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 803 END DO 804 END DO 805 806 IF ( ln_wd_diag ) THEN 807 WRITE(numout, *) 'zuwdmask: (i,j) = ', zuwdmask(iwdg,jwdg) 808 WRITE(numout, *) 'zwx: (i,j) = ', zwx(iwdg,jwdg) 809 WRITE(numout, *) 'e2u: (i,j) = ', e2u(iwdg,jwdg) 810 WRITE(numout, *) 'ua_e: (i,j) = ', ua_e(iwdg,jwdg) 811 WRITE(numout, *) 'un_e: (i,j) = ', un_e(iwdg,jwdg) 812 WRITE(numout, *) 'zhup2_e: (i,j) = ', zhup2_e(iwdg,jwdg) 813 WRITE(numout, *) 'zvwdmask: (i,j) = ', zvwdmask(iwdg,jwdg) 814 WRITE(numout, *) 'zwy: (i,j) = ', zwy(iwdg,jwdg) 815 END IF 816 817 END IF 818 710 819 ! Sum over sub-time-steps to compute advective velocities 711 820 za2 = wgtbtp2(jn) 712 821 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 713 822 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 714 ! 823 824 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_rwd_bc = True) 825 IF ( ln_rwd_bc ) THEN 826 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 827 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 828 829 IF ( ln_wd_diag ) THEN 830 WRITE(numout, *) 'za2, r1_e2u(i,j) = ', za2, r1_e2u(iwdg,jwdg) 831 WRITE(numout, *) 'un_adv: (i,j) = ', un_adv(iwdg,jwdg) 832 WRITE(numout, *) 'zuwdav2: (i,j) = ', zuwdav2(iwdg,jwdg) 833 WRITE(numout, *) 'zvwdav2: (i,j) = ', zvwdav2(iwdg,jwdg) 834 END IF 835 836 END IF 837 715 838 ! Set next sea level: 716 839 DO jj = 2, jpjm1 … … 770 893 DO ji = 2, jpim1 771 894 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 772 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &773 & MAX( zsshp2_e(ji,jj) + ht_ wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) &895 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 896 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj), zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 774 897 & > rn_wdmin1 + rn_wdmin2 775 898 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 776 899 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 777 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )900 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 778 901 779 902 IF(ll_tmp1) THEN … … 781 904 ELSE IF(ll_tmp2) THEN 782 905 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 783 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_ wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &906 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 784 907 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 785 908 ELSE … … 788 911 789 912 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 790 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &791 & MAX( zsshp2_e(ji,jj) + ht_ wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) &913 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 914 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj), zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 792 915 & > rn_wdmin1 + rn_wdmin2 793 916 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 794 917 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 795 & MAX( -ht_ wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )918 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 796 919 797 920 IF(ll_tmp1) THEN … … 799 922 ELSE IF(ll_tmp2) THEN 800 923 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 801 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_ wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &924 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 802 925 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 803 926 ELSE … … 886 1009 ! 887 1010 ! Add bottom stresses: 888 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 889 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 1011 !jth do implicitly instead 1012 ! zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 1013 ! zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 890 1014 ! 891 1015 ! Add top stresses: … … 933 1057 & + zv_frc(ji,jj) ) & 934 1058 & ) * ssvmask(ji,jj) 1059 1060 !jth implicit bottom friction: 1061 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * bfrua(ji,jj) * hur_e(ji,jj)) 1062 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * bfrva(ji,jj) * hvr_e(ji,jj)) 1063 935 1064 END DO 936 1065 END DO … … 941 1070 942 1071 IF( ln_wd ) THEN 943 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1)944 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1)1072 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 1073 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 945 1074 ELSE 946 1075 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) … … 964 1093 END DO 965 1094 ENDIF 966 ! 1095 1096 ! if ln_rwd: ua_e and va_e should not be masked ; they are used to determine the direction of flow into all cells 1097 1098 ! IF ( ln_rwd) THEN 1099 ! IF ( ln_wd_diag ) THEN 1100 ! WRITE(numout, *) 'ua_e: (i,j) = ', ua_e(iwdg,jwdg) 1101 ! WRITE(numout, *) 'va_e: (i,j) = ', va_e(iwdg,jwdg) 1102 ! END IF 1103 ! ua_e(:,:) = ua_e(:,:) * zuwdmask(:,:) 1104 ! va_e(:,:) = va_e(:,:) * zvwdmask(:,:) 1105 ! IF ( ln_wd_diag ) THEN 1106 ! WRITE(numout, *) 'ua_e: (i,j) = ', ua_e(iwdg,jwdg) 1107 ! WRITE(numout, *) 'va_e: (i,j) = ', va_e(iwdg,jwdg) 1108 ! END IF 1109 ! END IF 1110 1111 967 1112 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 968 1113 IF( ln_wd ) THEN 969 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1)970 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1)1114 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1115 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 971 1116 ELSE 972 1117 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) … … 1006 1151 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 1007 1152 ELSE ! Sum transports 1008 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 1009 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 1153 IF (.NOT.ln_rwd) THEN 1154 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 1155 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 1156 ELSE 1157 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 1158 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 1159 END IF 1010 1160 ENDIF 1011 1161 ! ! Sum sea level 1012 1162 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 1163 1013 1164 ! ! ==================== ! 1014 1165 END DO ! end loop ! … … 1033 1184 vb2_b(:,:) = zwy(:,:) 1034 1185 ENDIF 1186 1187 IF ( ln_wd_diag ) THEN 1188 WRITE(numout, *) 'ub2_b: (i,j) = ', ub2_b(iwdg,jwdg) 1189 WRITE(numout, *) 'r1_hu_n: (i,j) = ', r1_hu_n(iwdg,jwdg) 1190 WRITE(numout, *) 'zwx: (i,j) = ', zwx(iwdg,jwdg) 1191 WRITE(numout, *) 'un_adv: (i,j) = ', un_adv(iwdg,jwdg) 1192 END IF 1193 1035 1194 ! 1036 1195 ! Update barotropic trend: … … 1062 1221 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1063 1222 ENDIF 1064 ! 1223 1224 IF ( ln_wd_diag ) THEN 1225 WRITE(numout, *) 'ua_b: (i,j) A = ', ua_b(iwdg,jwdg) 1226 WRITE(numout, *) 'va_b: (i,j) B = ', va_b(iwdg,jwdg) 1227 END IF 1228 1229 ! temporary debugging code 1230 IF ( ln_wd_diag ) THEN 1231 WRITE(numout, *) 'ua: (i,j,k) B = ', ua(iwdg,jwdg,kwdg) 1232 WRITE(numout, *) 'ua_b: (i,j) B = ', ua_b(iwdg,jwdg) 1233 WRITE(numout, *) 'un: (i,j,k) = ', un(iwdg,jwdg,kwdg) 1234 WRITE(numout, *) 'un_b: (i,j) = ', un_b(iwdg,jwdg) 1235 WRITE(numout, *) 'un_adv: (i,j) = ', un_adv(iwdg,jwdg) 1236 WRITE(numout, *) 'va: (i,j,k) = ', va(iwdg,jwdg,kwdg) 1237 WRITE(numout, *) 'va_b: (i,j,k) = ', va_b(iwdg,jwdg) 1238 WRITE(numout, *) 'vn: (i,j,k) = ', vn(iwdg,jwdg,kwdg) 1239 WRITE(numout, *) 'vn_b: (i,j) = ', vn_b(iwdg,jwdg) 1240 WRITE(numout, *) 'vn_adv: (i,j) = ', vn_adv(iwdg,jwdg) 1241 END IF 1242 1243 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 1065 1244 DO jk = 1, jpkm1 1066 ! Correct velocities:1067 1245 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 1068 1246 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 1069 !1070 1247 END DO 1071 ! 1248 1249 IF ( ln_rwd .and. ln_rwd_bc) THEN 1250 DO jk = 1, jpkm1 1251 un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)) ) * umask(:,:,jk) 1252 vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)) ) * vmask(:,:,jk) 1253 END DO 1254 END IF 1255 1256 IF ( ln_wd_diag ) THEN 1257 WRITE(numout, *) 'ua: (i,j,k) = ', ua(iwdg,jwdg,kwdg) 1258 WRITE(numout, *) 'ua_b: (i,j,k) = ', ua_b(iwdg,jwdg) 1259 WRITE(numout, *) 'un: (i,j,k) = ', un(iwdg,jwdg,kwdg) 1260 WRITE(numout, *) 'va: (i,j,k) = ', va(iwdg,jwdg,kwdg) 1261 WRITE(numout, *) 'va_b: (i,j,k) = ', va_b(iwdg,jwdg) 1262 WRITE(numout, *) 'vn: (i,j,k) = ', vn(iwdg,jwdg,kwdg) 1263 END IF 1264 1072 1265 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 1073 1266 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current … … 1098 1291 CALL wrk_dealloc( jpi,jpj, zhf ) 1099 1292 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1293 IF( ln_rwd ) CALL wrk_dealloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 1100 1294 ! 1101 1295 IF ( ln_diatmb ) THEN -
branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r7646 r8865 1 1 MODULE wet_dry 2 3 !! includes updates to namelist namwad for diagnostic outputs of ROMS wetting and drying 4 2 5 !!============================================================================== 3 6 !! *** MODULE wet_dry *** … … 32 35 33 36 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_wd !: wetting and drying t-pt depths35 37 ! (can include negative depths) 38 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdramp, wdrampu, wdrampv ! for hpg limiting 36 39 37 40 LOGICAL, PUBLIC :: ln_wd !: Wetting/drying activation switch (T:on,F:off) 41 REAL(wp), PUBLIC :: rn_wdmin0 !: depth at which wetting/drying starts 42 LOGICAL, PUBLIC :: ln_rwd !: ROMS Wetting/drying activation switch (T:on,F:off) 38 43 REAL(wp), PUBLIC :: rn_wdmin1 !: minimum water depth on dried cells 39 REAL(wp), PUBLIC :: rn_wdmin2 !: toler rance of minimum water depth on dried cells44 REAL(wp), PUBLIC :: rn_wdmin2 !: tolerance of minimum water depth on dried cells 40 45 REAL(wp), PUBLIC :: rn_wdld !: land elevation below which wetting/drying 41 46 !: will be considered 42 47 INTEGER , PUBLIC :: nn_wdit !: maximum number of iteration for W/D limiter 48 LOGICAL, PUBLIC :: ln_rwd_bc !: ROMS scheme: True implies 3D velocities are set to the barotropic values at points 49 !: where the flow is from wet points on less than half the barotropic sub-steps 50 LOGICAL, PUBLIC :: ln_rwd_rmp !: use a ramp for the rwd flux limiter between 2 rn_wdmin1 and rn_wdmin1 (rather than a cut-off at rn_wdmin1) 51 REAL(wp), PUBLIC :: rn_ssh_ref !: height of z=0 with respect to the geoid; 52 REAL(wp), PUBLIC :: rn_ht_0 !: the height at which ht_0 = 0 53 54 LOGICAL, PUBLIC :: ln_wd_diag ! True implies wad diagnostic at chosen point are printed out 55 INTEGER , PUBLIC :: jn_wd_i, jn_wd_j, jn_wd_k ! indices at which diagnostic outputs are generated 43 56 44 57 PUBLIC wad_init ! initialisation routine called by step.F90 … … 58 71 !! ** input : - namwad namelist 59 72 !!---------------------------------------------------------------------- 60 NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 73 NAMELIST/namwad/ ln_wd, ln_rwd, rn_wdmin0, ln_rwd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit, ln_wd_diag, ln_rwd_bc, & 74 & rn_ht_0, jn_wd_i, jn_wd_j, jn_wd_k,ln_rwd_rmp 75 61 76 INTEGER :: ios ! Local integer output status for namelist read 62 77 INTEGER :: ierr ! Local integer status array allocation … … 78 93 WRITE(numout,*) '~~~~~~~~' 79 94 WRITE(numout,*) ' Namelist namwad' 80 WRITE(numout,*) ' Logical activation ln_wd = ', ln_wd 95 WRITE(numout,*) ' Logical for NOC wd scheme ln_wd = ', ln_wd 96 WRITE(numout,*) ' Logical for ROMS wd scheme ln_rwd = ', ln_rwd 97 WRITE(numout,*) ' Depth at which wet/drying starts rn_wdmin0 = ', rn_wdmin0 81 98 WRITE(numout,*) ' Minimum wet depth on dried cells rn_wdmin1 = ', rn_wdmin1 82 99 WRITE(numout,*) ' Tolerance of min wet depth rn_wdmin2 = ', rn_wdmin2 83 100 WRITE(numout,*) ' land elevation threshold rn_wdld = ', rn_wdld 84 101 WRITE(numout,*) ' Max iteration for W/D limiter nn_wdit = ', nn_wdit 102 WRITE(numout,*) ' Logical for WAD diagnostics ln_wd_diag = ', ln_wd_diag 103 WRITE(numout,*) ' T => baroclinic u,v = 0 at dry pts: ln_rwd_bc = ', ln_rwd_bc 104 WRITE(numout,*) ' use a ramp for rwd limiter: ln_rwd_rmp = ', ln_rwd_rmp 105 WRITE(numout,*) ' the height (z) at which ht_0 = 0:rn_ht_0 = ', rn_ht_0 106 WRITE(numout,*) ' i-index for diagnostic point jn_wd_i = ', jn_wd_i 107 WRITE(numout,*) ' j-index for diagnostic point jn_wd_j = ', jn_wd_j 108 WRITE(numout,*) ' k-index for diagnostic point jn_wd_k = ', jn_wd_k 85 109 ENDIF 86 110 ! 87 IF(ln_wd) THEN 88 ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj), STAT=ierr ) 111 IF(ln_wd .OR. ln_rwd) THEN 112 ALLOCATE( wdmask(jpi,jpj), STAT=ierr ) 113 ALLOCATE( wdramp(jpi,jpj), wdrampu(jpi,jpj), wdrampv(jpi,jpj), STAT=ierr ) 89 114 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 90 115 ENDIF … … 161 186 162 187 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 dry188 IF( ht_0(ji,jj)-rn_ssh_ref > zdepwd ) CYCLE ! and cells which are unlikely to dry 164 189 165 190 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & … … 168 193 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 169 194 170 zdep2 = ht_ wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1195 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 171 196 IF(zdep2 .le. 0._wp) THEN !add more safty, but not necessary 172 sshb1(ji,jj) = rn_wdmin1 - ht_ wd(ji,jj)197 sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 173 198 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 174 199 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp … … 180 205 END DO 181 206 207 208 ! slwa 209 ! HPG limiter from jholt 210 wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp) 211 ! write(6,*)'wdramp ',wdramp(10,10),wdramp(10,11) 212 !jth assume don't need a lbc_lnk here 213 DO jj = 1, jpjm1 214 DO ji = 1, jpim1 215 wdrampu(ji,jj) = min(wdramp(ji,jj),wdramp(ji+1,jj)) 216 wdrampv(ji,jj) = min(wdramp(ji,jj),wdramp(ji,jj+1)) 217 ENDDO 218 ENDDO 219 !wdrampu(:,:)=1.0_wp 220 !wdrampv(:,:)=1.0_wp 221 ! end HPG limiter 222 223 182 224 183 225 !! start limiter iterations … … 193 235 194 236 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 195 IF( ht_ wd(ji,jj) > zdepwd ) CYCLE237 IF( ht_0(ji,jj) > zdepwd ) CYCLE 196 238 197 239 ztmp = e1e2t(ji,jj) … … 203 245 204 246 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 205 zdep2 = ht_ wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)247 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 206 248 207 249 IF( zdep1 > zdep2 ) THEN … … 317 359 318 360 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 dry361 IF( ht_0(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry 320 362 321 363 zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + & … … 324 366 & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp) 325 367 326 zdep2 = ht_ wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1368 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 327 369 IF(zdep2 .le. 0._wp) THEN !add more safety, but not necessary 328 sshn_e(ji,jj) = rn_wdmin1 - ht_ wd(ji,jj)370 sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 329 371 IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp 330 372 IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp … … 348 390 349 391 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE 350 IF( ht_ wd(ji,jj) > zdepwd ) CYCLE392 IF( ht_0(ji,jj) > zdepwd ) CYCLE 351 393 352 394 ztmp = e1e2t(ji,jj) … … 358 400 359 401 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 360 zdep2 = ht_ wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)402 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 361 403 362 404 IF(zdep1 > zdep2) THEN
Note: See TracChangeset
for help on using the changeset viewer.
