Changeset 11771 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/domvvl.F90
- Timestamp:
- 2019-10-23T15:04:25+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/domvvl.F90
r10425 r11771 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 9 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 10 11 !!---------------------------------------------------------------------- 11 12 … … 13 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 14 15 !! dom_vvl_sf_nxt : Compute next vertical scale factors 15 !! dom_vvl_sf_ swp: Swap vertical scale factors and update the vertical grid16 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid 16 17 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 17 18 !! dom_vvl_rst : read/write restart file … … 37 38 PUBLIC dom_vvl_init ! called by domain.F90 38 39 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 PUBLIC dom_vvl_sf_ swp! called by step.F9040 PUBLIC dom_vvl_sf_update ! called by step.F90 40 41 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 42 … … 93 94 94 95 95 SUBROUTINE dom_vvl_init 96 SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 96 97 !!---------------------------------------------------------------------- 97 98 !! *** ROUTINE dom_vvl_init *** … … 104 105 !! 105 106 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 106 !! - Regrid: e3 (u/v)_n107 !! e3 (u/v)_b108 !! e3w _n109 !! e3 (u/v)w_b110 !! e3 (u/v)w_n111 !! gdept _n, gdepw_n and gde3w_n107 !! - Regrid: e3[u/v](:,:,:,Kmm) 108 !! e3[u/v](:,:,:,Kmm) 109 !! e3w(:,:,:,Kmm) 110 !! e3[u/v]w_b 111 !! e3[u/v]w_n 112 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 112 113 !! - h(t/u/v)_0 113 114 !! - frq_rst_e3t and frq_rst_hdv … … 115 116 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 117 !!---------------------------------------------------------------------- 118 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 119 ! 117 120 INTEGER :: ji, jj, jk 118 121 INTEGER :: ii0, ii1, ij0, ij1 … … 130 133 ! 131 134 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 132 CALL dom_vvl_rst( nit000, 'READ' )133 e3t _a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all135 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 136 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 134 137 ! 135 138 ! !== Set of all other vertical scale factors ==! (now and before) 136 139 ! ! Horizontal interpolation of e3t 137 CALL dom_vvl_interpol( e3t _b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U138 CALL dom_vvl_interpol( e3t _n(:,:,:), e3u_n(:,:,:), 'U' )139 CALL dom_vvl_interpol( e3t _b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V140 CALL dom_vvl_interpol( e3t _n(:,:,:), e3v_n(:,:,:), 'V' )141 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F140 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U 141 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 142 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V 143 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 144 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F 142 145 ! ! Vertical interpolation of e3t,u,v 143 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W144 CALL dom_vvl_interpol( e3t _b(:,:,:), e3w_b (:,:,:), 'W' )145 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW146 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )147 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW148 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )146 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W 147 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' ) 148 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW 149 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 150 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW 151 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 149 152 150 153 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 151 e3t _a(:,:,:) = e3t_n(:,:,:)152 e3u _a(:,:,:) = e3u_n(:,:,:)153 e3v _a(:,:,:) = e3v_n(:,:,:)154 e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 155 e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 156 e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 154 157 ! 155 158 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 156 gdept _n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration)157 gdepw _n(:,:,1) = 0.0_wp158 gde3w _n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg159 gdept _b(:,:,1) = 0.5_wp * e3w_b(:,:,1)160 gdepw _b(:,:,1) = 0.0_wp159 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) ! reference to the ocean surface (used for MLD and light penetration) 160 gdepw(:,:,1,Kmm) = 0.0_wp 161 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg 162 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 163 gdepw(:,:,1,Kbb) = 0.0_wp 161 164 DO jk = 2, jpk ! vertical sum 162 165 DO jj = 1,jpj … … 165 168 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 166 169 ! ! 0.5 where jk = mikt 167 !!gm ??????? BUG ? gdept _n as well as gde3w_ndoes not include the thickness of ISF ??170 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 168 171 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 169 gdepw _n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)170 gdept _n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &171 & + (1-zcoef) * ( gdept _n(ji,jj,jk-1) + e3w_n(ji,jj,jk))172 gde3w _n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)173 gdepw _b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)174 gdept _b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) &175 & + (1-zcoef) * ( gdept _b(ji,jj,jk-1) + e3w_b(ji,jj,jk))172 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 173 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 174 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 175 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 176 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 177 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 178 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 176 179 END DO 177 180 END DO … … 179 182 ! 180 183 ! !== thickness of the water column !! (ocean portion only) 181 ht _n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....182 hu _b(:,:) = e3u_b(:,:,1) * umask(:,:,1)183 hu _n(:,:) = e3u_n(:,:,1) * umask(:,:,1)184 hv _b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)185 hv _n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)184 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 185 hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 186 hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 187 hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 188 hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 186 189 DO jk = 2, jpkm1 187 ht _n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)188 hu _b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)189 hu _n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)190 hv _b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)191 hv _n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)190 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 191 hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 192 hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 193 hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 194 hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 192 195 END DO 193 196 ! 194 197 ! !== inverse of water column thickness ==! (u- and v- points) 195 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF196 r1_hu _n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )197 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )198 r1_hv _n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )198 r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 199 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 200 r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 201 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 199 202 200 203 ! !== z_tilde coordinate case ==! (Restoring frequencies) … … 266 269 267 270 268 SUBROUTINE dom_vvl_sf_nxt( kt, kcall )271 SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 269 272 !!---------------------------------------------------------------------- 270 273 !! *** ROUTINE dom_vvl_sf_nxt *** … … 288 291 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 289 292 !!---------------------------------------------------------------------- 290 INTEGER, INTENT( in ) :: kt ! time step 291 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 293 INTEGER, INTENT( in ) :: kt ! time step 294 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step 295 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 292 296 ! 293 297 INTEGER :: ji, jj, jk ! dummy loop indices … … 321 325 ! ! --------------------------------------------- ! 322 326 ! 323 z_scale(:,:) = ( ssh a(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )327 z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 324 328 DO jk = 1, jpkm1 325 ! formally this is the same as e3t _a= e3t_0*(1+ssha/ht_0)326 e3t _a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)329 ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 330 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 327 331 END DO 328 332 ! … … 337 341 zht(:,:) = 0._wp 338 342 DO jk = 1, jpkm1 339 zhdiv(:,:) = zhdiv(:,:) + e3t _n(:,:,jk) * hdivn(:,:,jk)340 zht (:,:) = zht (:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)343 zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 344 zht (:,:) = zht (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 341 345 END DO 342 346 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 348 352 DO jk = 1, jpkm1 349 353 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 350 & * ( hdiv_lf(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )354 & * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 351 355 END DO 352 356 ENDIF … … 361 365 IF( ln_vvl_ztilde ) THEN ! z_tilde case 362 366 DO jk = 1, jpkm1 363 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )367 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 364 368 END DO 365 369 ELSE ! layer case 366 370 DO jk = 1, jpkm1 367 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)371 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 368 372 END DO 369 373 ENDIF … … 476 480 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 477 481 END DO 478 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh n(:,:) + 1. - ssmask(:,:) )482 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 479 483 DO jk = 1, jpkm1 480 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t _n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)484 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 481 485 END DO 482 486 … … 486 490 ! ! ---baroclinic part--------- ! 487 491 DO jk = 1, jpkm1 488 e3t _a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)492 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 489 493 END DO 490 494 ENDIF … … 501 505 zht(:,:) = 0.0_wp 502 506 DO jk = 1, jpkm1 503 zht(:,:) = zht(:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)504 END DO 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh n(:,:) - zht(:,:) ) )507 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 508 END DO 509 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 506 510 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t _n))) =', z_tmax511 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 508 512 ! 509 513 zht(:,:) = 0.0_wp 510 514 DO jk = 1, jpkm1 511 zht(:,:) = zht(:,:) + e3t _a(:,:,jk) * tmask(:,:,jk)512 END DO 513 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh a(:,:) - zht(:,:) ) )515 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 516 END DO 517 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 514 518 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 515 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t _a))) =', z_tmax519 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 516 520 ! 517 521 zht(:,:) = 0.0_wp 518 522 DO jk = 1, jpkm1 519 zht(:,:) = zht(:,:) + e3t _b(:,:,jk) * tmask(:,:,jk)520 END DO 521 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh b(:,:) - zht(:,:) ) )523 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 524 END DO 525 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 522 526 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 523 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t _b))) =', z_tmax524 ! 525 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh b(:,:) ) )527 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 528 ! 529 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) ) 526 530 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 527 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh b))) =', z_tmax528 ! 529 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh n(:,:) ) )531 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 532 ! 533 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) ) 530 534 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 531 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh n))) =', z_tmax532 ! 533 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh a(:,:) ) )535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 536 ! 537 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) ) 534 538 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh a))) =', z_tmax539 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 536 540 END IF 537 541 … … 540 544 ! *********************************** ! 541 545 542 CALL dom_vvl_interpol( e3t _a(:,:,:), e3u_a(:,:,:), 'U' )543 CALL dom_vvl_interpol( e3t _a(:,:,:), e3v_a(:,:,:), 'V' )546 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 547 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 544 548 545 549 ! *********************************** ! … … 547 551 ! *********************************** ! 548 552 549 hu _a(:,:) = e3u_a(:,:,1) * umask(:,:,1)550 hv _a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)553 hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 554 hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 551 555 DO jk = 2, jpkm1 552 hu _a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)553 hv _a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)556 hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 557 hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 554 558 END DO 555 559 ! ! Inverse of the local depth 556 560 !!gm BUG ? don't understand the use of umask_i here ..... 557 r1_hu _a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )558 r1_hv _a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )561 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 562 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 559 563 ! 560 564 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 563 567 564 568 565 SUBROUTINE dom_vvl_sf_ swp( kt)566 !!---------------------------------------------------------------------- 567 !! *** ROUTINE dom_vvl_sf_ swp***569 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 570 !!---------------------------------------------------------------------- 571 !! *** ROUTINE dom_vvl_sf_update *** 568 572 !! 569 !! ** Purpose : compute time filter and swap of scale factors573 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 570 574 !! compute all depths and related variables for next time step 571 575 !! write outputs and restart file 572 576 !! 573 !! ** Method : - swap of e3t with trick for volume/tracer conservation 577 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 574 578 !! - reconstruct scale factor at other grid points (interpolate) 575 579 !! - recompute depths and water height fields 576 580 !! 577 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_nready for next time step581 !! ** Action : - tilde_e3t_(b/n) ready for next time step 578 582 !! - Recompute: 579 583 !! e3(u/v)_b 580 !! e3w _n584 !! e3w(:,:,:,Kmm) 581 585 !! e3(u/v)w_b 582 586 !! e3(u/v)w_n 583 !! gdept _n, gdepw_n and gde3w_n587 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 584 588 !! h(u/v) and h(u/v)r 585 589 !! … … 587 591 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 588 592 !!---------------------------------------------------------------------- 589 INTEGER, INTENT( in ) :: kt ! time step 593 INTEGER, INTENT( in ) :: kt ! time step 594 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 590 595 ! 591 596 INTEGER :: ji, jj, jk ! dummy loop indices … … 595 600 IF( ln_linssh ) RETURN ! No calculation in linear free surface 596 601 ! 597 IF( ln_timing ) CALL timing_start('dom_vvl_sf_ swp')602 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 598 603 ! 599 604 IF( kt == nit000 ) THEN 600 605 IF(lwp) WRITE(numout,*) 601 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_ swp : - time filter and swap of scale factors'602 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'606 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 607 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 603 608 ENDIF 604 609 ! … … 615 620 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 616 621 ENDIF 617 gdept _b(:,:,:) = gdept_n(:,:,:)618 gdepw _b(:,:,:) = gdepw_n(:,:,:)619 620 e3t _n(:,:,:) = e3t_a(:,:,:)621 e3u _n(:,:,:) = e3u_a(:,:,:)622 e3v _n(:,:,:) = e3v_a(:,:,:)622 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 623 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 624 625 e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa) 626 e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa) 627 e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa) 623 628 624 629 ! Compute all missing vertical scale factor and depths … … 626 631 ! Horizontal scale factor interpolations 627 632 ! -------------------------------------- 628 ! - ML - e3u _b and e3v_b are allready computed in dynnxt629 ! - JC - hu _b, hv_b, hur_b, hvr_b also633 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 634 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 630 635 631 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' )636 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 632 637 633 638 ! Vertical scale factor interpolations 634 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n(:,:,:), 'W' )635 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' )636 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' )637 CALL dom_vvl_interpol( e3t _b(:,:,:), e3w_b(:,:,:), 'W' )638 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )639 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )639 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 640 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 641 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 642 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 643 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 644 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 640 645 641 646 ! t- and w- points depth (set the isf depth as it is in the initial step) 642 gdept _n(:,:,1) = 0.5_wp * e3w_n(:,:,1)643 gdepw _n(:,:,1) = 0.0_wp644 gde3w _n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)647 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 648 gdepw(:,:,1,Kmm) = 0.0_wp 649 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 645 650 DO jk = 2, jpk 646 651 DO jj = 1,jpj … … 649 654 ! 1 for jk = mikt 650 655 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 651 gdepw _n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)652 gdept _n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) &653 & + (1-zcoef) * ( gdept _n(ji,jj,jk-1) + e3w_n(ji,jj,jk) )654 gde3w _n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)656 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 657 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 658 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 659 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 655 660 END DO 656 661 END DO … … 659 664 ! Local depth and Inverse of the local depth of the water 660 665 ! ------------------------------------------------------- 661 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 662 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 663 ! 664 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 666 ! 667 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 665 668 DO jk = 2, jpkm1 666 ht _n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)669 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 667 670 END DO 668 671 669 672 ! write restart file 670 673 ! ================== 671 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )672 ! 673 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_ swp')674 ! 675 END SUBROUTINE dom_vvl_sf_ swp674 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 675 ! 676 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 677 ! 678 END SUBROUTINE dom_vvl_sf_update 676 679 677 680 … … 783 786 784 787 785 SUBROUTINE dom_vvl_rst( kt, cdrw )788 SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 786 789 !!--------------------------------------------------------------------- 787 790 !! *** ROUTINE dom_vvl_rst *** … … 795 798 !! they are set to 0. 796 799 !!---------------------------------------------------------------------- 797 INTEGER , INTENT(in) :: kt ! ocean time-step 798 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 800 INTEGER , INTENT(in) :: kt ! ocean time-step 801 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 802 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 799 803 ! 800 804 INTEGER :: ji, jj, jk … … 806 810 IF( ln_rstart ) THEN !* Read the restart file 807 811 CALL rst_read_open ! open the restart file if necessary 808 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh n, ldxios = lrxios )812 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios ) 809 813 ! 810 814 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) … … 817 821 ! ! --------- ! 818 822 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 819 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t _b(:,:,:), ldxios = lrxios )820 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t _n(:,:,:), ldxios = lrxios )823 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 824 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 821 825 ! needed to restart if land processor not computed 822 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t _b and e3t_nfound in restart files'826 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 823 827 WHERE ( tmask(:,:,:) == 0.0_wp ) 824 e3t _n(:,:,:) = e3t_0(:,:,:)825 e3t _b(:,:,:) = e3t_0(:,:,:)828 e3t(:,:,:,Kmm) = e3t_0(:,:,:) 829 e3t(:,:,:,Kbb) = e3t_0(:,:,:) 826 830 END WHERE 827 831 IF( neuler == 0 ) THEN 828 e3t _b(:,:,:) = e3t_n(:,:,:)832 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 829 833 ENDIF 830 834 ELSE IF( id1 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart files'835 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 832 836 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 833 837 IF(lwp) write(numout,*) 'neuler is forced to 0' 834 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t _b(:,:,:), ldxios = lrxios )835 e3t _n(:,:,:) = e3t_b(:,:,:)838 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 839 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 836 840 neuler = 0 837 841 ELSE IF( id2 > 0 ) THEN 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _bnot found in restart files'842 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 839 843 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 840 844 IF(lwp) write(numout,*) 'neuler is forced to 0' 841 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t _n(:,:,:), ldxios = lrxios )842 e3t _b(:,:,:) = e3t_n(:,:,:)845 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 846 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 843 847 neuler = 0 844 848 ELSE 845 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart file'849 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 846 850 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 847 851 IF(lwp) write(numout,*) 'neuler is forced to 0' 848 852 DO jk = 1, jpk 849 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &853 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 850 854 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 851 855 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 852 856 END DO 853 e3t _b(:,:,:) = e3t_n(:,:,:)857 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 854 858 neuler = 0 855 859 ENDIF … … 888 892 IF( cn_cfg == 'wad' ) THEN 889 893 ! Wetting and drying test case 890 CALL usr_def_istate( gdept _b, tmask, tsb, ub, vb, sshb)891 ts n (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones892 ssh n (:,:) = sshb(:,:)893 u n (:,:,:) = ub (:,:,:)894 v n (:,:,:) = vb (:,:,:)894 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 895 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 896 ssh (:,:,Kmm) = ssh(:,:,Kbb) 897 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 898 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 895 899 ELSE 896 900 ! if not test case 897 ssh n(:,:) = -ssh_ref898 ssh b(:,:) = -ssh_ref901 ssh(:,:,Kmm) = -ssh_ref 902 ssh(:,:,Kbb) = -ssh_ref 899 903 900 904 DO jj = 1, jpj 901 905 DO ji = 1, jpi 902 906 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 903 904 sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 905 sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 906 ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 907 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 908 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 907 909 ENDIF 908 910 ENDDO … … 912 914 ! Adjust vertical metrics for all wad 913 915 DO jk = 1, jpk 914 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &916 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 915 917 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 916 918 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 917 919 END DO 918 e3t _b(:,:,:) = e3t_n(:,:,:)920 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 919 921 920 922 DO ji = 1, jpi … … 928 930 ELSE 929 931 ! 930 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t _b and e3t_n931 CALL usr_def_istate( gdept_0, tmask, ts b, ub, vb, sshb )932 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm) 933 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 932 934 ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 933 935 ! 934 936 DO jk=1,jpk 935 e3t _b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &936 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &937 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t_b != 0 on land points937 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 938 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 939 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t_b != 0 on land points 938 940 END DO 939 e3t_n(:,:,:) = e3t_b(:,:,:) 940 sshn(:,:) = sshb(:,:) ! needed later for gde3w 941 !!$ e3t_n(:,:,:)=e3t_0(:,:,:) 942 !!$ e3t_b(:,:,:)=e3t_0(:,:,:) 941 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 942 ssh(:,: ,Kmm) = ssh(:,: ,Kbb) ! needed later for gde3w 943 943 ! 944 944 END IF ! end of ll_wd edits … … 958 958 ! ! all cases ! 959 959 ! ! --------- ! 960 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t _b(:,:,:), ldxios = lwxios )961 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t _n(:,:,:), ldxios = lwxios )960 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 961 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 962 962 ! ! ----------------------- ! 963 963 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
Note: See TracChangeset
for help on using the changeset viewer.
