- Timestamp:
- 2020-04-03T18:54:55+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90
r12679 r12680 9 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 10 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 11 !! 4.x ! 2020-02 (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate11 !! 4.x ! 2020-02 (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate 12 12 !!---------------------------------------------------------------------- 13 13 14 14 !!---------------------------------------------------------------------- 15 !! dom_qe_init 16 !! dom_qe_sf_nxt 17 !! dom_qe_sf_update 15 !! dom_qe_init : define initial vertical scale factors, depths and column thickness 16 !! dom_qe_sf_nxt : Compute next vertical scale factors 17 !! dom_qe_sf_update: Swap vertical scale factors and update the vertical grid 18 18 !! dom_qe_interpol : Interpolate vertical scale factors from one grid point to another 19 !! dom_qe_r3c 20 !! dom_qe_rst: read/write restart file21 !! dom_qe_ctl 19 !! dom_qe_r3c : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points 20 !! qe_rst_read : read/write restart file 21 !! dom_qe_ctl : Check the vvl options 22 22 !!---------------------------------------------------------------------- 23 USE oce 24 USE phycst 25 USE dom_oce 23 USE oce ! ocean dynamics and tracers 24 USE phycst ! physical constant 25 USE dom_oce ! ocean space and time domain 26 26 USE dynadv , ONLY : ln_dynadv_vec 27 USE isf_oce 28 USE sbc_oce 29 USE wet_dry 30 USE usrdef_istate 31 USE restart 27 USE isf_oce ! iceshelf cavities 28 USE sbc_oce ! ocean surface boundary condition 29 USE wet_dry ! wetting and drying 30 USE usrdef_istate ! user defined initial state (wad only) 31 USE restart ! ocean restart 32 32 ! 33 USE in_out_manager 34 USE iom 35 USE lib_mpp 36 USE lbclnk 37 USE timing 33 USE in_out_manager ! I/O manager 34 USE iom ! I/O manager library 35 USE lib_mpp ! distributed memory computing library 36 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 37 USE timing ! Timing 38 38 39 39 IMPLICIT NONE … … 42 42 PUBLIC dom_qe_init ! called by domain.F90 43 43 PUBLIC dom_qe_zgr ! called by isfcpl.F90 44 PUBLIC dom_qe_sf_nxt ! called by steplf.F9045 PUBLIC dom_qe_sf_update ! called by steplf.F9044 !!st PUBLIC dom_qe_sf_nxt ! called by steplf.F90 45 !!st PUBLIC dom_qe_sf_update ! called by steplf.F90 46 46 PUBLIC dom_h_nxt ! called by steplf.F90 47 PUBLIC dom_h_update ! called by steplf.F90 47 48 PUBLIC dom_qe_r3c ! called by steplf.F90 48 49 … … 102 103 ! 103 104 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 104 CALL dom_qe_rst( nit000, Kbb, Kmm, 'READ')105 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all105 CALL qe_rst_read( nit000, Kbb, Kmm ) 106 !!st e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 106 107 ! 107 108 CALL dom_qe_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 108 109 ! 109 IF(lwxios) THEN ! define variables in restart file when writing with XIOS110 CALL iom_set_rstw_var_active('e3t_b')111 CALL iom_set_rstw_var_active('e3t_n')112 ENDIF110 ! IF(lwxios) THEN ! define variables in restart file when writing with XIOS 111 ! CALL iom_set_rstw_var_active('e3t_b') 112 ! CALL iom_set_rstw_var_active('e3t_n') 113 ! ENDIF 113 114 ! 114 115 END SUBROUTINE dom_qe_init … … 147 148 CALL dom_qe_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) ) 148 149 ! 149 DO jk = 1, jpkm1 ! Horizontal interpolation of e3t 150 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) * tmask(:,:,jk) ) ! Kbb time level 151 e3u(:,:,jk,Kbb) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) * umask(:,:,jk) ) 152 e3v(:,:,jk,Kbb) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) * vmask(:,:,jk) ) 153 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) ) ! Kmm time level 154 e3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) * umask(:,:,jk) ) 155 e3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) * vmask(:,:,jk) ) 156 e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) ) 157 END DO 158 ! 159 DO jk = 1, jpk ! Vertical interpolation of e3t,u,v 160 ! ! The ratio does not have to be masked at w-level 161 e3w (:,:,jk,Kbb) = e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) ) ! Kbb time level 162 e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) ) 163 e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) ) 164 e3w (:,:,jk,Kmm) = e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) ) ! Kmm time level 165 e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) ) 166 e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) ) 167 END DO 168 ! 169 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 170 e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 171 e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 172 e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 150 ! !!st 151 ! DO jk = 1, jpkm1 ! Horizontal interpolation of e3t 152 ! e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) * tmask(:,:,jk) ) ! Kbb time level 153 ! e3u(:,:,jk,Kbb) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) * umask(:,:,jk) ) 154 ! e3v(:,:,jk,Kbb) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) * vmask(:,:,jk) ) 155 ! e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) ) ! Kmm time level 156 ! e3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) * umask(:,:,jk) ) 157 ! e3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) * vmask(:,:,jk) ) 158 ! e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) ) 159 ! END DO 160 ! ! 161 ! DO jk = 1, jpk ! Vertical interpolation of e3t,u,v 162 ! ! ! The ratio does not have to be masked at w-level 163 ! e3w (:,:,jk,Kbb) = e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) ) ! Kbb time level 164 ! e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) ) 165 ! e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) ) 166 ! e3w (:,:,jk,Kmm) = e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) ) ! Kmm time level 167 ! e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) ) 168 ! e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) ) 169 ! END DO 170 ! ! 171 ! ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 172 ! e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 173 ! e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 174 ! e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 175 !!st end 173 176 ! 174 177 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) … … 221 224 END SUBROUTINE dom_qe_zgr 222 225 223 224 SUBROUTINE dom_qe_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 226 ! !!st 227 ! SUBROUTINE dom_qe_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 228 ! !!---------------------------------------------------------------------- 229 ! !! *** ROUTINE dom_qe_sf_nxt *** 230 ! !! 231 ! !! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt, 232 ! !! tranxt and dynspg routines 233 ! !! 234 ! !! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness. 235 ! !! 236 ! !! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case 237 ! !! - tilde_e3t_a: after increment of vertical scale factor 238 ! !! in z_tilde case 239 ! !! - e3(t/u/v)_a 240 ! !! 241 ! !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 242 ! !!---------------------------------------------------------------------- 243 ! INTEGER, INTENT( in ) :: kt ! time step 244 ! INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step 245 ! INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 246 ! ! 247 ! INTEGER :: ji, jj, jk ! dummy loop indices 248 ! INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 249 ! REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars 250 ! LOGICAL :: ll_do_bclinic ! local logical 251 ! REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv 252 ! !!---------------------------------------------------------------------- 253 ! ! 254 ! IF( ln_linssh ) RETURN ! No calculation in linear free surface 255 ! ! 256 ! IF( ln_timing ) CALL timing_start('dom_qe_sf_nxt') 257 ! ! 258 ! IF( kt == nit000 ) THEN 259 ! IF(lwp) WRITE(numout,*) 260 ! IF(lwp) WRITE(numout,*) 'dom_qe_sf_nxt : compute after scale factors' 261 ! IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 262 ! ENDIF 263 ! 264 ! 265 ! ! ******************************* ! 266 ! ! After acale factors at t-points ! 267 ! ! ******************************* ! 268 ! ! ! --------------------------------------------- ! 269 ! ! ! z_star coordinate and barotropic z-tilde part ! 270 ! ! ! --------------------------------------------- ! 271 ! ! 272 ! ! 273 ! ! *********************************** ! 274 ! ! After scale factors at u- v- points ! 275 ! ! *********************************** ! 276 ! ! 277 ! DO jk = 1, jpkm1 278 ! e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) * tmask(:,:,jk) ) 279 ! e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) * umask(:,:,jk) ) 280 ! e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) * vmask(:,:,jk) ) 281 ! END DO 282 ! ! 283 ! ! *********************************** ! 284 ! ! After depths at u- v points ! 285 ! ! *********************************** ! 286 ! hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) ) 287 ! hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) ) 288 ! ! ! Inverse of the local depth 289 ! r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 290 ! r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 291 ! ! 292 ! IF( ln_timing ) CALL timing_stop('dom_qe_sf_nxt') 293 ! ! 294 ! END SUBROUTINE dom_qe_sf_nxt 295 !!st end 296 297 SUBROUTINE dom_h_nxt( kt, Kbb, Kmm, Kaa, kcall ) 225 298 !!---------------------------------------------------------------------- 226 299 !! *** ROUTINE dom_qe_sf_nxt *** 227 300 !! 228 !! ** Purpose : - compute the after scale factorsused in tra_zdf, dynnxt,301 !! ** Purpose : - compute the after water heigh used in tra_zdf, dynnxt, 229 302 !! tranxt and dynspg routines 230 303 !! 231 !! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness. 232 !! 233 !! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case 234 !! - tilde_e3t_a: after increment of vertical scale factor 235 !! in z_tilde case 236 !! - e3(t/u/v)_a 237 !! 238 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 304 !! ** Method : - z_star case: Proportionnaly to the water column thickness. 305 !! 306 !! ** Action : - h(u/v) update wrt ssh/h(u/v)_0 307 !! 239 308 !!---------------------------------------------------------------------- 240 309 INTEGER, INTENT( in ) :: kt ! time step … … 242 311 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 243 312 ! 244 INTEGER :: ji, jj, jk ! dummy loop indices245 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers246 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars247 LOGICAL :: ll_do_bclinic ! local logical248 REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv249 313 !!---------------------------------------------------------------------- 250 314 ! 251 315 IF( ln_linssh ) RETURN ! No calculation in linear free surface 252 316 ! 253 IF( ln_timing ) CALL timing_start('dom_ qe_sf_nxt')317 IF( ln_timing ) CALL timing_start('dom_h_nxt') 254 318 ! 255 319 IF( kt == nit000 ) THEN 256 320 IF(lwp) WRITE(numout,*) 257 IF(lwp) WRITE(numout,*) 'dom_ qe_sf_nxt : compute after scale factors'321 IF(lwp) WRITE(numout,*) 'dom_h_nxt : compute after scale factors' 258 322 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 259 323 ENDIF 260 261 262 ! ******************************* !263 ! After acale factors at t-points !264 ! ******************************* !265 ! ! --------------------------------------------- !266 ! ! z_star coordinate and barotropic z-tilde part !267 ! ! --------------------------------------------- !268 !269 !270 ! *********************************** !271 ! After scale factors at u- v- points !272 ! *********************************** !273 !274 DO jk = 1, jpkm1275 e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) * tmask(:,:,jk) )276 e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) * umask(:,:,jk) )277 e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) * vmask(:,:,jk) )278 END DO279 324 ! 280 325 ! *********************************** ! … … 287 332 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 288 333 ! 289 IF( ln_timing ) CALL timing_stop('dom_qe_sf_nxt')290 !291 END SUBROUTINE dom_qe_sf_nxt292 293 294 SUBROUTINE dom_h_nxt( kt, Kbb, Kmm, Kaa, kcall )295 !!----------------------------------------------------------------------296 !! *** ROUTINE dom_qe_sf_nxt ***297 !!298 !! ** Purpose : - compute the after water heigh used in tra_zdf, dynnxt,299 !! tranxt and dynspg routines300 !!301 !! ** Method : - z_star case: Proportionnaly to the water column thickness.302 !!303 !! ** Action : - h(u/v) update wrt ssh/h(u/v)_0304 !!305 !!----------------------------------------------------------------------306 INTEGER, INTENT( in ) :: kt ! time step307 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step308 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence309 !310 !!----------------------------------------------------------------------311 !312 IF( ln_linssh ) RETURN ! No calculation in linear free surface313 !314 IF( ln_timing ) CALL timing_start('dom_h_nxt')315 !316 IF( kt == nit000 ) THEN317 IF(lwp) WRITE(numout,*)318 IF(lwp) WRITE(numout,*) 'dom_h_nxt : compute after scale factors'319 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'320 ENDIF321 !322 ! *********************************** !323 ! After depths at u- v points !324 ! *********************************** !325 hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) )326 hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) )327 ! ! Inverse of the local depth328 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )329 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )330 !331 334 IF( ln_timing ) CALL timing_stop('dom_h_nxt') 332 335 ! 333 336 END SUBROUTINE dom_h_nxt 334 337 335 336 SUBROUTINE dom_qe_sf_update( kt, Kbb, Kmm, Kaa ) 338 ! !!st 339 ! SUBROUTINE dom_qe_sf_update( kt, Kbb, Kmm, Kaa ) 340 ! !!---------------------------------------------------------------------- 341 ! !! *** ROUTINE dom_qe_sf_update *** 342 ! !! 343 ! !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 344 ! !! compute all depths and related variables for next time step 345 ! !! write outputs and restart file 346 ! !! 347 ! !! ** Method : - reconstruct scale factor at other grid points (interpolate) 348 ! !! - recompute depths and water height fields 349 ! !! 350 ! !! ** Action : - Recompute: 351 ! !! e3(u/v)_b 352 ! !! e3w(:,:,:,Kmm) 353 ! !! e3(u/v)w_b 354 ! !! e3(u/v)w_n 355 ! !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 356 ! !! h(u/v) and h(u/v)r 357 ! !! 358 ! !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 359 ! !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 360 ! !!---------------------------------------------------------------------- 361 ! INTEGER, INTENT( in ) :: kt ! time step 362 ! INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 363 ! ! 364 ! INTEGER :: ji, jj, jk ! dummy loop indices 365 ! REAL(wp) :: zcoef ! local scalar 366 ! !!---------------------------------------------------------------------- 367 ! ! 368 ! IF( ln_linssh ) RETURN ! No calculation in linear free surface 369 ! ! 370 ! IF( ln_timing ) CALL timing_start('dom_qe_sf_update') 371 ! ! 372 ! IF( kt == nit000 ) THEN 373 ! IF(lwp) WRITE(numout,*) 374 ! IF(lwp) WRITE(numout,*) 'dom_qe_sf_update : - interpolate scale factors and compute depths for next time step' 375 ! IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 376 ! ENDIF 377 ! ! 378 ! ! Compute all missing vertical scale factor and depths 379 ! ! ==================================================== 380 ! ! Horizontal scale factor interpolations 381 ! ! -------------------------------------- 382 ! ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 383 ! ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 384 ! 385 ! 386 ! ! Scale factor computation 387 ! DO jk = 1, jpk ! Horizontal interpolation 388 ! e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) ) ! Kmm time level 389 ! ! ! Vertical interpolation 390 ! ! ! The ratio does not have to be masked at w-level 391 ! e3w (:,:,jk,Kmm) = e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) ) ! Kmm time level 392 ! e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) ) 393 ! e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) ) 394 ! e3w (:,:,jk,Kbb) = e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) ) ! Kbb time level 395 ! e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) ) 396 ! e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) ) 397 ! END DO 398 ! 399 ! 400 ! IF( ln_isf ) THEN !** IceShelF cavities 401 ! ! ! to be created depending of the new names in isf 402 ! ! ! it should be something like that : (with h_isf = thickness of iceshelf) 403 ! ! ! in fact currently, h_isf(:,:) is called : risfdep(:,:) 404 ! !!gm - depth idea 0 : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask ! 405 ! gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) ) 406 ! gdepw(:,:,1,Kmm) = 0._wp ! Initialized to zero one for all 407 ! gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg 408 ! DO jk = 2, jpk 409 ! gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) ) & 410 ! + MAX( 0._wp , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) ) 411 ! gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) ) & 412 ! + MAX( 0._wp , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) ) 413 ! gde3w(:,:,jk) = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm) 414 ! gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) ) & 415 ! + MAX( 0._wp , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) ) 416 ! gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) ) & 417 ! + MAX( 0._wp , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) ) 418 ! END DO 419 ! ! 420 ! ELSE !** No cavities (all depth rescaled, even inside topography: no mask) 421 ! ! 422 ! !!gm idea 0 : just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask ! 423 ! DO jk = 1, jpk 424 ! gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) ) 425 ! gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) ) 426 ! gde3w(:,:,jk) = gdept (:,:,jk,Kmm) - ssh(:,:,Kmm) 427 ! gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) ) 428 ! gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) ) 429 ! END DO 430 ! ! 431 ! ENDIF 432 ! 433 ! ! Local depth and Inverse of the local depth of the water 434 ! ! ------------------------------------------------------- 435 ! ! 436 ! ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm) 437 ! 438 ! ! write restart file 439 ! ! ================== 440 ! IF( lrst_oce ) CALL dom_qe_rst( kt, Kbb, Kmm, 'WRITE' ) 441 ! ! 442 ! IF( ln_timing ) CALL timing_stop('dom_qe_sf_update') 443 ! ! 444 ! END SUBROUTINE dom_qe_sf_update 445 !!st end 446 447 SUBROUTINE dom_h_update( kt, Kbb, Kmm, Kaa ) 337 448 !!---------------------------------------------------------------------- 338 449 !! *** ROUTINE dom_qe_sf_update *** … … 346 457 !! 347 458 !! ** Action : - Recompute: 348 !! e3(u/v)_b349 !! e3w(:,:,:,Kmm)350 !! e3(u/v)w_b351 !! e3(u/v)w_n352 459 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 353 460 !! h(u/v) and h(u/v)r … … 377 484 ! Horizontal scale factor interpolations 378 485 ! -------------------------------------- 379 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt380 486 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 381 382 383 ! Scale factor computation384 DO jk = 1, jpk ! Horizontal interpolation385 e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) ) ! Kmm time level386 ! ! Vertical interpolation387 ! ! The ratio does not have to be masked at w-level388 e3w (:,:,jk,Kmm) = e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) ) ! Kmm time level389 e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )390 e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )391 e3w (:,:,jk,Kbb) = e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) ) ! Kbb time level392 e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )393 e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )394 END DO395 396 487 397 488 IF( ln_isf ) THEN !** IceShelF cavities … … 399 490 ! ! it should be something like that : (with h_isf = thickness of iceshelf) 400 491 ! ! in fact currently, h_isf(:,:) is called : risfdep(:,:) 401 !!gm - depth idea 0 : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !492 !!gm - depth idea 0 : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask ! 402 493 gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) ) 403 494 gdepw(:,:,1,Kmm) = 0._wp ! Initialized to zero one for all … … 417 508 ELSE !** No cavities (all depth rescaled, even inside topography: no mask) 418 509 ! 419 !!gm idea 0 : just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !510 !!gm idea 0 : just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask ! 420 511 DO jk = 1, jpk 421 512 gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) ) … … 435 526 ! write restart file 436 527 ! ================== 437 IF( lrst_oce ) CALL dom_qe_rst( kt, Kbb, Kmm, 'WRITE' )438 !439 528 IF( ln_timing ) CALL timing_stop('dom_qe_sf_update') 440 529 ! 441 END SUBROUTINE dom_ qe_sf_update530 END SUBROUTINE dom_h_update 442 531 443 532 … … 507 596 508 597 509 SUBROUTINE dom_qe_rst( kt, Kbb, Kmm, cdrw)598 SUBROUTINE qe_rst_read( kt, Kbb, Kmm ) 510 599 !!--------------------------------------------------------------------- 511 !! *** ROUTINE dom_qe_rst***512 !! 513 !! ** Purpose : Read or write VVL filein restart file600 !! *** ROUTINE qe_rst_read *** 601 !! 602 !! ** Purpose : Read ssh in restart file 514 603 !! 515 604 !! ** Method : use of IOM library 516 !! if the restart does not contain vertical scale factors, 517 !! they are set to the _0 values 518 !! if the restart does not contain vertical scale factors increments (z_tilde), 519 !! they are set to 0. 605 !! if the restart does not contain ssh, 606 !! it is set to the _0 values. 520 607 !!---------------------------------------------------------------------- 521 608 INTEGER , INTENT(in) :: kt ! ocean time-step 522 609 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 523 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag524 610 ! 525 611 INTEGER :: ji, jj, jk … … 527 613 !!---------------------------------------------------------------------- 528 614 ! 529 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise530 ! ! ===============531 615 IF( ln_rstart ) THEN !* Read the restart file 532 616 CALL rst_read_open ! open the restart file if necessary 533 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios )534 617 ! 535 id1 = iom_varid( numror, ' e3t_b', ldstop = .FALSE. )536 id2 = iom_varid( numror, ' e3t_n', ldstop = .FALSE. )618 id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. ) 619 id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. ) 537 620 ! 538 621 ! ! --------- ! … … 541 624 ! 542 625 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 543 CALL iom_get( numror, jpdom_autoglo, ' e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios)544 CALL iom_get( numror, jpdom_autoglo, ' e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios)626 CALL iom_get( numror, jpdom_autoglo, 'sshb' , ssh(:,:,Kbb), ldxios = lrxios ) 627 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios ) 545 628 ! needed to restart if land processor not computed 546 IF(lwp) write(numout,*) ' dom_qe_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'547 WHERE ( tmask(:,:,:) == 0.0_wp )548 e3t(:,:,:,Kmm) = e3t_0(:,:,:)549 e3t(:,:,:,Kbb) = e3t_0(:,:,:)629 IF(lwp) write(numout,*) 'qe_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files' 630 WHERE ( ssmask(:,:) == 0.0_wp ) !!gm/st ==> sm should not be necessary on ssh when it was required on e3 631 ssh(:,:,Kmm) = 0._wp 632 ssh(:,:,Kbb) = 0._wp 550 633 END WHERE 551 634 IF( neuler == 0 ) THEN 552 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)635 ssh(:,:,Kbb) = ssh(:,:,Kmm) 553 636 ENDIF 554 637 ELSE IF( id1 > 0 ) THEN 555 IF(lwp) write(numout,*) ' dom_qe_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'556 IF(lwp) write(numout,*) ' e3t_n set equal to e3t_b.'638 IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart files' 639 IF(lwp) write(numout,*) 'sshn set equal to sshb.' 557 640 IF(lwp) write(numout,*) 'neuler is forced to 0' 558 CALL iom_get( numror, jpdom_autoglo, ' e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )559 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)641 CALL iom_get( numror, jpdom_autoglo, 'sshb', ssh(:,:,Kbb), ldxios = lrxios ) 642 ssh(:,:,Kmm) = ssh(:,:,Kbb) 560 643 neuler = 0 561 644 ELSE IF( id2 > 0 ) THEN 562 IF(lwp) write(numout,*) ' dom_qe_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'563 IF(lwp) write(numout,*) ' e3t_b set equal to e3t_n.'645 IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kbb) not found in restart files' 646 IF(lwp) write(numout,*) 'sshb set equal to sshn.' 564 647 IF(lwp) write(numout,*) 'neuler is forced to 0' 565 CALL iom_get( numror, jpdom_autoglo, ' e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )566 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)648 CALL iom_get( numror, jpdom_autoglo, 'sshn', ssh(:,:,Kmm), ldxios = lrxios ) 649 ssh(:,:,Kbb) = ssh(:,:,Kmm) 567 650 neuler = 0 568 651 ELSE 569 IF(lwp) write(numout,*) ' dom_qe_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'570 IF(lwp) write(numout,*) ' Compute scale factor from sshn'652 IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart file' 653 IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero' 571 654 IF(lwp) write(numout,*) 'neuler is forced to 0' 572 DO jk = 1, jpk 573 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 574 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 575 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 576 END DO 577 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 655 ssh(:,:,:) = 0._wp 578 656 neuler = 0 579 657 ENDIF … … 583 661 IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential 584 662 ! 585 IF( cn_cfg == 'wad' ) THEN 586 ! Wetting and drying test case 663 IF( cn_cfg == 'wad' ) THEN ! Wetting and drying test case 587 664 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 588 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 589 ssh (:,:,Kmm) = ssh(:,:,Kbb) 590 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 591 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 592 ELSE 593 ! if not test case 665 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 666 ssh(:,: ,Kmm) = ssh(:,: ,Kbb) 667 uu (:,:,: ,Kmm) = uu (:,:,: ,Kbb) 668 vv (:,:,: ,Kmm) = vv (:,:,: ,Kbb) 669 ELSE ! if not test case 594 670 ssh(:,:,Kmm) = -ssh_ref 595 671 ssh(:,:,Kbb) = -ssh_ref 596 672 ! 597 673 DO_2D_11_11 598 674 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth … … 601 677 ENDIF 602 678 END_2D 603 ENDIF !If test case else 604 605 ! Adjust vertical metrics for all wad 606 DO jk = 1, jpk 607 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) ) 608 END DO 609 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 679 ENDIF 610 680 611 681 DO ji = 1, jpi 612 682 DO jj = 1, jpj 613 683 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 614 CALL ctl_stop( ' dom_qe_rst: ht_0 must be positive at potentially wet points' )684 CALL ctl_stop( 'qe_rst_read: ht_0 must be positive at potentially wet points' ) 615 685 ENDIF 616 686 END DO … … 623 693 ! CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 624 694 ! ! 625 ! DO jk=1,jpk 626 ! e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 627 ! & / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 628 ! END DO 629 ! e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 630 ssh(:,:,Kmm)=0._wp 631 e3t(:,:,:,Kmm)=e3t_0(:,:,:) 632 e3t(:,:,:,Kbb)=e3t_0(:,:,:) 695 ssh(:,:,:) = 0._wp 633 696 ! 634 697 ENDIF ! end of ll_wd edits 635 698 ! 636 699 ENDIF 637 ! 638 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 639 ! ! =================== 640 IF(lwp) WRITE(numout,*) '---- dom_qe_rst ----' 641 IF( lwxios ) CALL iom_swap( cwxios_context ) 642 ! ! --------- ! 643 ! ! all cases ! 644 ! ! --------- ! 645 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 646 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 647 ! 648 IF( lwxios ) CALL iom_swap( cxios_context ) 649 ENDIF 650 ! 651 END SUBROUTINE dom_qe_rst 700 ! 701 END SUBROUTINE qe_rst_read 652 702 653 703
Note: See TracChangeset
for help on using the changeset viewer.
