Changeset 11412 for NEMO/branches
- Timestamp:
- 2019-08-06T17:57:38+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/VORTEX/MY_SRC/domvvl.F90
r10572 r11412 93 93 94 94 95 SUBROUTINE dom_vvl_init 95 SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 96 96 !!---------------------------------------------------------------------- 97 97 !! *** ROUTINE dom_vvl_init *** … … 104 104 !! 105 105 !! ** 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_n106 !! - Regrid: e3[u/v](:,:,:,Kmm) 107 !! e3[u/v](:,:,:,Kmm) 108 !! e3w(:,:,:,Kmm) 109 !! e3[u/v]w_b 110 !! e3[u/v]w_n 111 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 112 112 !! - h(t/u/v)_0 113 113 !! - frq_rst_e3t and frq_rst_hdv … … 115 115 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 116 !!---------------------------------------------------------------------- 117 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 118 ! 117 119 INTEGER :: ji, jj, jk 118 120 INTEGER :: ii0, ii1, ij0, ij1 … … 130 132 ! 131 133 ! ! 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 all134 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 135 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 134 136 ! 135 137 ! !== Set of all other vertical scale factors ==! (now and before) 136 138 ! ! 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 F139 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U 140 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 141 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V 142 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 143 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F 142 144 ! ! 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' )145 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W 146 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' ) 147 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW 148 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 149 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW 150 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 149 151 150 152 ! 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(:,:,:)153 e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 154 e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 155 e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 154 156 ! 155 157 ! !== 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_wp158 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) ! reference to the ocean surface (used for MLD and light penetration) 159 gdepw(:,:,1,Kmm) = 0.0_wp 160 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg 161 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 162 gdepw(:,:,1,Kbb) = 0.0_wp 161 163 DO jk = 2, jpk ! vertical sum 162 164 DO jj = 1,jpj … … 165 167 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 166 168 ! ! 0.5 where jk = mikt 167 !!gm ??????? BUG ? gdept _n as well as gde3w_ndoes not include the thickness of ISF ??169 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 168 170 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))171 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 172 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 173 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 174 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 175 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 176 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 177 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 176 178 END DO 177 179 END DO … … 179 181 ! 180 182 ! !== 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)183 ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 184 hu_b(:,:) = e3u(:,:,1,Kbb) * umask(:,:,1) 185 hu_n(:,:) = e3u(:,:,1,Kmm) * umask(:,:,1) 186 hv_b(:,:) = e3v(:,:,1,Kbb) * vmask(:,:,1) 187 hv_n(:,:) = e3v(:,:,1,Kmm) * vmask(:,:,1) 186 188 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)189 ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 190 hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 191 hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 192 hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 193 hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 192 194 END DO 193 195 ! … … 266 268 267 269 268 SUBROUTINE dom_vvl_sf_nxt( kt, kcall )270 SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 269 271 !!---------------------------------------------------------------------- 270 272 !! *** ROUTINE dom_vvl_sf_nxt *** … … 288 290 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 289 291 !!---------------------------------------------------------------------- 290 INTEGER, INTENT( in ) :: kt ! time step 291 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 292 INTEGER, INTENT( in ) :: kt ! time step 293 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step 294 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 292 295 ! 293 296 INTEGER :: ji, jj, jk ! dummy loop indices … … 321 324 ! ! --------------------------------------------- ! 322 325 ! 323 z_scale(:,:) = ( ssh a(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )326 z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 324 327 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)328 ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 329 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 327 330 END DO 328 331 ! … … 337 340 zht(:,:) = 0._wp 338 341 DO jk = 1, jpkm1 339 zhdiv(:,:) = zhdiv(:,:) + e3t _n(:,:,jk) * hdivn(:,:,jk)340 zht (:,:) = zht (:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)342 zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 343 zht (:,:) = zht (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 341 344 END DO 342 345 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 348 351 DO jk = 1, jpkm1 349 352 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 350 & * ( hdiv_lf(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )353 & * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 351 354 END DO 352 355 ENDIF … … 361 364 IF( ln_vvl_ztilde ) THEN ! z_tilde case 362 365 DO jk = 1, jpkm1 363 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )366 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 364 367 END DO 365 368 ELSE ! layer case 366 369 DO jk = 1, jpkm1 367 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)370 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 368 371 END DO 369 372 ENDIF … … 476 479 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 477 480 END DO 478 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh n(:,:) + 1. - ssmask(:,:) )481 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 479 482 DO jk = 1, jpkm1 480 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t _n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)483 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 481 484 END DO 482 485 … … 486 489 ! ! ---baroclinic part--------- ! 487 490 DO jk = 1, jpkm1 488 e3t _a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)491 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 489 492 END DO 490 493 ENDIF … … 501 504 zht(:,:) = 0.0_wp 502 505 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(:,:) ) )506 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 507 END DO 508 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 506 509 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_tmax510 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 508 511 ! 509 512 zht(:,:) = 0.0_wp 510 513 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(:,:) ) )514 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 515 END DO 516 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 514 517 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_tmax518 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 516 519 ! 517 520 zht(:,:) = 0.0_wp 518 521 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(:,:) ) )522 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 523 END DO 524 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 522 525 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(:,:) ) )526 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 527 ! 528 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) ) 526 529 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(:,:) ) )530 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 531 ! 532 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) ) 530 533 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(:,:) ) )534 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 535 ! 536 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) ) 534 537 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh a))) =', z_tmax538 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 536 539 END IF 537 540 … … 540 543 ! *********************************** ! 541 544 542 CALL dom_vvl_interpol( e3t _a(:,:,:), e3u_a(:,:,:), 'U' )543 CALL dom_vvl_interpol( e3t _a(:,:,:), e3v_a(:,:,:), 'V' )545 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 546 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 544 547 545 548 ! *********************************** ! … … 547 550 ! *********************************** ! 548 551 549 hu_a(:,:) = e3u _a(:,:,1) * umask(:,:,1)550 hv_a(:,:) = e3v _a(:,:,1) * vmask(:,:,1)552 hu_a(:,:) = e3u(:,:,1,Kaa) * umask(:,:,1) 553 hv_a(:,:) = e3v(:,:,1,Kaa) * vmask(:,:,1) 551 554 DO jk = 2, jpkm1 552 hu_a(:,:) = hu_a(:,:) + e3u _a(:,:,jk) * umask(:,:,jk)553 hv_a(:,:) = hv_a(:,:) + e3v _a(:,:,jk) * vmask(:,:,jk)555 hu_a(:,:) = hu_a(:,:) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 556 hv_a(:,:) = hv_a(:,:) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 554 557 END DO 555 558 ! ! Inverse of the local depth … … 563 566 564 567 565 SUBROUTINE dom_vvl_sf_swp( kt )568 SUBROUTINE dom_vvl_sf_swp( kt, Kbb, Kmm, Kaa ) 566 569 !!---------------------------------------------------------------------- 567 570 !! *** ROUTINE dom_vvl_sf_swp *** … … 578 581 !! - Recompute: 579 582 !! e3(u/v)_b 580 !! e3w _n583 !! e3w(:,:,:,Kmm) 581 584 !! e3(u/v)w_b 582 585 !! e3(u/v)w_n 583 !! gdept _n, gdepw_n and gde3w_n586 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 584 587 !! h(u/v) and h(u/v)r 585 588 !! … … 587 590 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 588 591 !!---------------------------------------------------------------------- 589 INTEGER, INTENT( in ) :: kt ! time step 592 INTEGER, INTENT( in ) :: kt ! time step 593 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 590 594 ! 591 595 INTEGER :: ji, jj, jk ! dummy loop indices … … 615 619 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 616 620 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(:,:,:)621 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 622 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 623 624 e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa) 625 e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa) 626 e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa) 623 627 624 628 ! Compute all missing vertical scale factor and depths … … 626 630 ! Horizontal scale factor interpolations 627 631 ! -------------------------------------- 628 ! - ML - e3u _b and e3v_bare allready computed in dynnxt632 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are allready computed in dynnxt 629 633 ! - JC - hu_b, hv_b, hur_b, hvr_b also 630 634 631 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' )635 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 632 636 633 637 ! 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' )638 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 639 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 640 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 641 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 642 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 643 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 640 644 641 645 ! 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(:,:)646 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 647 gdepw(:,:,1,Kmm) = 0.0_wp 648 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 645 649 DO jk = 2, jpk 646 650 DO jj = 1,jpj … … 649 653 ! 1 for jk = mikt 650 654 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)655 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 656 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 657 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 658 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 655 659 END DO 656 660 END DO … … 662 666 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 663 667 ! 664 ht_n(:,:) = e3t _n(:,:,1) * tmask(:,:,1)668 ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 665 669 DO jk = 2, jpkm1 666 ht_n(:,:) = ht_n(:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)670 ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 667 671 END DO 668 672 669 673 ! write restart file 670 674 ! ================== 671 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )675 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 672 676 ! 673 677 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_swp') … … 783 787 784 788 785 SUBROUTINE dom_vvl_rst( kt, cdrw )789 SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 786 790 !!--------------------------------------------------------------------- 787 791 !! *** ROUTINE dom_vvl_rst *** … … 795 799 !! they are set to 0. 796 800 !!---------------------------------------------------------------------- 797 INTEGER , INTENT(in) :: kt ! ocean time-step 798 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 801 INTEGER , INTENT(in) :: kt ! ocean time-step 802 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 803 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 799 804 ! 800 805 INTEGER :: ji, jj, jk … … 806 811 IF( ln_rstart ) THEN !* Read the restart file 807 812 CALL rst_read_open ! open the restart file if necessary 808 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh n, ldxios = lrxios )813 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios ) 809 814 ! 810 815 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) … … 817 822 ! ! --------- ! 818 823 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 )824 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 825 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 821 826 ! needed to restart if land processor not computed 822 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t _b and e3t_nfound in restart files'827 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 823 828 WHERE ( tmask(:,:,:) == 0.0_wp ) 824 e3t _n(:,:,:) = e3t_0(:,:,:)825 e3t _b(:,:,:) = e3t_0(:,:,:)829 e3t(:,:,:,Kmm) = e3t_0(:,:,:) 830 e3t(:,:,:,Kbb) = e3t_0(:,:,:) 826 831 END WHERE 827 832 IF( neuler == 0 ) THEN 828 e3t _b(:,:,:) = e3t_n(:,:,:)833 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 829 834 ENDIF 830 835 ELSE IF( id1 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart files'836 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 832 837 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 833 838 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(:,:,:)839 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 840 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 836 841 neuler = 0 837 842 ELSE IF( id2 > 0 ) THEN 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _bnot found in restart files'843 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 839 844 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 840 845 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(:,:,:)846 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 847 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 843 848 neuler = 0 844 849 ELSE 845 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart file'850 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 846 851 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 847 852 IF(lwp) write(numout,*) 'neuler is forced to 0' 848 853 DO jk = 1, jpk 849 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &854 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 850 855 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 851 856 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 852 857 END DO 853 e3t _b(:,:,:) = e3t_n(:,:,:)858 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 854 859 neuler = 0 855 860 ENDIF … … 888 893 IF( cn_cfg == 'wad' ) THEN 889 894 ! 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 (:,:,:)895 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 896 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 897 ssh (:,:,Kmm) = ssh(:,:,Kbb) 898 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 899 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 895 900 ELSE 896 901 ! if not test case 897 ssh n(:,:) = -ssh_ref898 ssh b(:,:) = -ssh_ref902 ssh(:,:,Kmm) = -ssh_ref 903 ssh(:,:,Kbb) = -ssh_ref 899 904 900 905 DO jj = 1, jpj 901 906 DO ji = 1, jpi 902 907 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) ) 908 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 909 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 907 910 ENDIF 908 911 ENDDO … … 912 915 ! Adjust vertical metrics for all wad 913 916 DO jk = 1, jpk 914 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &917 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 915 918 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 916 919 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 917 920 END DO 918 e3t _b(:,:,:) = e3t_n(:,:,:)921 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 919 922 920 923 DO ji = 1, jpi … … 928 931 ELSE 929 932 ! 930 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n 931 CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb ) 933 ! ! 934 ! usr_def_istate called here only to get ssh, that is needed to initialize e3t_b and e3t_n 935 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 932 936 ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 933 937 ! 934 938 DO jk=1,jpk 935 e3t _b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) )&936 & 937 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t_b!= 0 on land points939 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 940 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 941 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t(Kbb) != 0 on land points 938 942 END DO 939 e3t_n(:,:,:) = e3t_b(:,:,:) 940 sshn(:,:) = sshb(:,:) ! needed later for gde3w 941 !!$ e3t_n(:,:,:)=e3t_0(:,:,:) 942 !!$ e3t_b(:,:,:)=e3t_0(:,:,:) 943 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 944 ssh(:,: ,Kmm) = ssh(:,: ,Kbb) ! needed later for gde3w 943 945 ! 944 946 END IF ! end of ll_wd edits … … 958 960 ! ! all cases ! 959 961 ! ! --------- ! 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 )962 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 963 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 962 964 ! ! ----------------------- ! 963 965 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
Note: See TracChangeset
for help on using the changeset viewer.
