Changeset 10001 for NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM
- Timestamp:
- 2018-07-26T09:50:51+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM
- Files:
-
- 5 edited
-
dom_oce.F90 (modified) (4 diffs)
-
domain.F90 (modified) (12 diffs)
-
dommsk.F90 (modified) (1 diff)
-
domvvl.F90 (modified) (16 diffs)
-
iscplrst.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90
r9939 r10001 30 30 !! ---------------------------- 31 31 ! !!* Namelist namdom : time & space domain * 32 LOGICAL , PUBLIC :: ln_RK3 !: 2rd order Runge-Kutta scheme 33 LOGICAL , PUBLIC :: ln_MLF !: Modified Leap-Frog scheme (Leclair & Madec OM2010) 34 REAL(wp), PUBLIC :: rn_atfp !: asselin time filter parameter 35 REAL(wp), PUBLIC :: rn_Dt !: time step for the dynamics and tracer (for both RK3 & MLF) 32 36 LOGICAL , PUBLIC :: ln_linssh !: =T linear free surface ==>> model level are fixed in time 33 37 LOGICAL , PUBLIC :: ln_meshmask !: =T create a mesh-mask file (mesh_mask.nc) 34 38 REAL(wp), PUBLIC :: rn_isfhmin !: threshold to discriminate grounded ice to floating ice 35 REAL(wp), PUBLIC :: rn_dt !: time step for the dynamics and tracer36 REAL(wp), PUBLIC :: rn_atfp !: asselin time filter parameter37 39 LOGICAL , PUBLIC :: ln_1st_euler !: =0 start with forward time step or not (=1) 38 40 LOGICAL , PUBLIC :: ln_iscpl !: coupling with ice sheet … … 143 145 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hu_b , hu_n , hu_a !: u-depth [m] 144 146 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0 , hv_b , hv_n , hv_a !: v-depth [m] 145 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_b , r1_hu_n , r1_hu_a !: inverse of u-depth [1/m] 146 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_b , r1_hv_n , r1_hv_a !: inverse of v-depth [1/m] 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hf_0 !: v-depth [m] 148 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ht_0 !: inverse of u-depth [1/m] 149 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a !: inverse of u-depth [1/m] 150 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a !: inverse of v-depth [1/m] 151 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hf_0 !: inverse of v-depth [1/m] 147 152 148 153 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) 149 154 INTEGER, PUBLIC :: nlb10 !: shallowest W level Bellow ~10m (nla10 + 1) 150 155 151 !! 1D reference 152 !! = -----------------====------156 !! 1D reference vertical coordinate 157 !! ==------------------------------ 153 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m) 154 159 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m) … … 170 175 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep !: Iceshelf draft (ISF) 171 176 172 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask 177 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts 173 178 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 174 179 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask !: land/ocean mask at WT-, WU- and WV-pts … … 265 270 & e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , STAT=ierr(5) ) 266 271 ! 267 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,&268 & hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,&269 & ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,&270 & hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6) )271 !272 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , hf_0(jpi,jpj) , & 273 & hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) , & 274 & ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) , & 275 & hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , & 276 & r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj), r1_hv_0(jpi,jpj) , r1_hf_0(jpi,jpj) , STAT=ierr(6) ) 272 277 ! 273 278 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(7) ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90
r9939 r10001 78 78 CHARACTER (len=*), INTENT(IN) :: cdstr ! model: NEMO or SAS. Determines core restart variables 79 79 INTEGER , DIMENSION(jpi,jpj) :: ik_top , ik_bot ! top and bottom ocean level 80 REAL(wp), DIMENSION(jpi,jpj) :: z1_hu_0, z1_hv_081 80 !!---------------------------------------------------------------------- 82 81 ! … … 138 137 CALL dom_msk( ik_top, ik_bot ) ! Masks 139 138 IF( ln_closea ) CALL dom_clo ! ln_closea=T : closed seas included in the simulation 140 139 ! ! Read in masks to define closed seas and lakes 141 140 ! 142 141 DO jj = 1, jpj ! depth of the iceshelves 143 142 DO ji = 1, jpi 144 143 ik = mikt(ji,jj) 145 risfdep(ji,jj) = gdepw_0(ji,jj,ik) 144 risfdep(ji,jj) = gdepw_0(ji,jj,ik) !!gm RENAME it as h_isf(:,:) better no? 146 145 END DO 147 146 END DO 148 147 ! 149 ht_0(:,:) = 0._wp ! Reference oceanthickness150 hu_0(:,:) = 0._wp 148 ht_0(:,:) = 0._wp ! Reference water column thickness (wet point only, i.e. without ISF thickness 149 hu_0(:,:) = 0._wp ! ------------ and zero over land ) 151 150 hv_0(:,:) = 0._wp 151 hf_0(:,:) = 0._wp 152 152 DO jk = 1, jpk 153 153 ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 154 154 hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 155 155 hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 156 hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * fmask(:,:,jk) 156 157 END DO 157 158 ! 159 r1_ht_0(:,:) = ssmask(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) ! =0 on lands 160 r1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) ) ! not in ISF areas 161 r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 162 r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp - ssfmask(:,:) ) 163 ! 158 164 ! !== time varying part of coordinate system ==! 159 165 ! 160 166 IF( ln_linssh ) THEN != Fix in time : set to the reference one for all 161 ! 162 ! before ! now ! after ! 163 gdept_b = gdept_0 ; gdept_n = gdept_0 ! --- ! depth of grid-points 164 gdepw_b = gdepw_0 ; gdepw_n = gdepw_0 ! --- ! 165 gde3w_n = gde3w_0 ! --- ! 166 ! 167 e3t_b = e3t_0 ; e3t_n = e3t_0 ; e3t_a = e3t_0 ! scale factors 168 e3u_b = e3u_0 ; e3u_n = e3u_0 ; e3u_a = e3u_0 ! 169 e3v_b = e3v_0 ; e3v_n = e3v_0 ; e3v_a = e3v_0 ! 170 e3f_n = e3f_0 ! --- ! 171 e3w_b = e3w_0 ; e3w_n = e3w_0 ! --- ! 172 e3uw_b = e3uw_0 ; e3uw_n = e3uw_0 ! --- ! 173 e3vw_b = e3vw_0 ; e3vw_n = e3vw_0 ! --- ! 174 ! 175 z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 176 z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 177 ! 178 ! before ! now ! after ! 179 ht_n = ht_0 ! ! water column thickness 180 hu_b = hu_0 ; hu_n = hu_0 ; hu_a = hu_0 ! 181 hv_b = hv_0 ; hv_n = hv_0 ; hv_a = hv_0 ! 182 r1_hu_b = z1_hu_0 ; r1_hu_n = z1_hu_0 ; r1_hu_a = z1_hu_0 ! inverse of water column thickness 183 r1_hv_b = z1_hv_0 ; r1_hv_n = z1_hv_0 ; r1_hv_a = z1_hv_0 ! 167 ! 168 ! before ! now ! after ! 169 gdept_b = gdept_0 ; gdept_n = gdept_0 ! --- ! depth 170 gdepw_b = gdepw_0 ; gdepw_n = gdepw_0 ! --- ! of 171 gde3w_n = gde3w_0 ! --- ! grid-points 172 ! ! ! ! 173 e3t_b = e3t_0 ; e3t_n = e3t_0 ; e3t_a = e3t_0 ! scale 174 e3u_b = e3u_0 ; e3u_n = e3u_0 ; e3u_a = e3u_0 ! factors 175 e3v_b = e3v_0 ; e3v_n = e3v_0 ; e3v_a = e3v_0 ! 176 e3f_n = e3f_0 ! --- ! 177 e3w_b = e3w_0 ; e3w_n = e3w_0 ! --- ! 178 e3uw_b = e3uw_0 ; e3uw_n = e3uw_0 ! --- ! 179 e3vw_b = e3vw_0 ; e3vw_n = e3vw_0 ! --- ! 180 ! ! 181 ht_n = ht_0 ! ! water column 182 hu_b = hu_0 ; hu_n = hu_0 ; hu_a = hu_0 ! thickness 183 hv_b = hv_0 ; hv_n = hv_0 ; hv_a = hv_0 ! 184 r1_hu_b = r1_hu_0 ; r1_hu_n = r1_hu_0 ; r1_hu_a = r1_hu_0 ! 1 / water column 185 r1_hv_b = r1_hv_0 ; r1_hv_n = r1_hv_0 ; r1_hv_a = r1_hv_0 ! thickness 184 186 ! 185 187 ! … … 192 194 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 193 195 ! 194 IF( ln_meshmask .AND. .NOT.ln_iscpl )CALL dom_wri ! Create a domain file196 IF( ln_meshmask .AND. .NOT.ln_iscpl ) CALL dom_wri ! Create a domain file 195 197 IF( ln_meshmask .AND. ln_iscpl .AND. .NOT.ln_rstart ) CALL dom_wri ! Create a domain file 196 198 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 197 199 ! 198 IF( ln_write_cfg ) CALL cfg_write 200 IF( ln_write_cfg ) CALL cfg_write ! create the configuration file 199 201 ! 200 202 IF(lwp) THEN … … 204 206 WRITE(numout,*) 205 207 ENDIF 206 !208 ! 207 209 END SUBROUTINE dom_init 208 210 … … 290 292 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, & 291 293 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl , & 292 & nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate , &294 & nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , ln_istate , & 293 295 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, ln_1st_euler, & 294 296 & ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 295 NAMELIST/namdom/ ln_ linssh, rn_isfhmin, rn_Dt, rn_atfp, ln_crs, ln_meshmask297 NAMELIST/namdom/ ln_RK3, ln_MLF, rn_Dt, rn_atfp, ln_linssh, rn_isfhmin, ln_crs, ln_meshmask 296 298 #if defined key_netcdf4 297 299 NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip … … 305 307 ENDIF 306 308 ! 309 REWIND( numnam_ref ) ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep) 310 READ ( numnam_ref, namdom, IOSTAT = ios, ERR = 903) 311 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp ) 312 REWIND( numnam_cfg ) ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep) 313 READ ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 314 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 315 IF(lwm) WRITE( numond, namdom ) 316 ! 317 IF(lwp) THEN 318 WRITE(numout,*) 319 WRITE(numout,*) ' Namelist : namdom --- space & time domain' 320 WRITE(numout,*) ' 3rd order Runge-Kutta scheme ln_RK3 = ', ln_RK3 321 WRITE(numout,*) ' Modified Leap-Frog scheme ln_MLF = ', ln_MLF 322 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 323 WRITE(numout,*) ' ocean time step (RK3 & MLF) rn_Dt = ', rn_Dt 324 WRITE(numout,*) ' linear free surface (=T) ln_linssh = ', ln_linssh 325 WRITE(numout,*) ' create mesh/mask file ln_meshmask = ', ln_meshmask 326 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' [m]' 327 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs 328 WRITE(numout,*) 329 ENDIF 330 ! 331 IF ( ( ln_RK3 .AND. .NOT.ln_MLF ) .OR. & 332 & ( .NOT.ln_RK3 .AND. ln_MLF ) ) THEN 333 IF( ln_RK3 .AND. lwp ) WRITE(numout,*)' ==>>> a RK3 time-stepping scheme is used' 334 IF( ln_MLF .AND. lwp ) WRITE(numout,*)' ==>>> a MLF time-stepping scheme is used' 335 ELSE 336 CALL ctl_stop( 'dom_nam: select one and only one time stepping scheme') 337 ENDIF 307 338 ! 308 339 REWIND( numnam_ref ) ! Namelist namrun in reference namelist : Parameters of the run … … 315 346 ! 316 347 IF(lwp) THEN ! control print 348 WRITE(numout,*) 317 349 WRITE(numout,*) ' Namelist : namrun --- run parameters' 318 350 WRITE(numout,*) ' Assimilation cycle nn_no = ', nn_no … … 330 362 WRITE(numout,*) ' initial time of day in hhmm nn_time0 = ', nn_time0 331 363 WRITE(numout,*) ' leap year calendar (0/1) nn_leapy = ', nn_leapy 332 WRITE(numout,*) ' initial state output nn_istate = ', nn_istate364 WRITE(numout,*) ' output the initial state ln_istate = ', ln_istate 333 365 IF( ln_rst_list ) THEN 334 366 WRITE(numout,*) ' list of restart dump times nn_stocklist =', nn_stocklist … … 357 389 ndate0 = nn_date0 358 390 nleapy = nn_leapy 359 ninist = nn_istate360 391 nstock = nn_stock 361 392 nstocklist = nn_stocklist 362 393 nwrite = nn_write 363 IF( ln_rstart ) THEN ! restart : set 1st time-step scheme (LF or forward) 364 l_1st_euler = ln_1st_euler 365 ELSE ! start from rest : always an Euler scheme for the 1st time-step 366 IF(lwp) WRITE(numout,*) 367 IF(lwp) WRITE(numout,*)' ==>>> Start from rest (ln_rstart=F)' 368 IF(lwp) WRITE(numout,*)' an Euler initial time step is used ' 369 l_1st_euler = .TRUE. 394 395 IF( ln_MLF ) THEN ! Leap-Frog only 396 IF( ln_rstart ) THEN ! restart : set 1st time-step scheme (LF or forward) 397 l_1st_euler = ln_1st_euler 398 ELSE ! start from rest : always an Euler scheme for the 1st time-step 399 IF(lwp) WRITE(numout,*) 400 IF(lwp) WRITE(numout,*)' ==>>> Start from rest (ln_rstart=F)' 401 IF(lwp) WRITE(numout,*)' an Euler initial time step is used ' 402 l_1st_euler = .TRUE. 403 ENDIF 370 404 ENDIF 371 405 ! ! control of output frequency … … 399 433 ENDIF 400 434 #endif 401 402 REWIND( numnam_ref ) ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep) 403 READ ( numnam_ref, namdom, IOSTAT = ios, ERR = 903) 404 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp ) 405 REWIND( numnam_cfg ) ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep) 406 READ ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 407 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 408 IF(lwm) WRITE( numond, namdom ) 409 ! 410 IF(lwp) THEN 411 WRITE(numout,*) 412 WRITE(numout,*) ' Namelist : namdom --- space & time domain' 413 WRITE(numout,*) ' linear free surface (=T) ln_linssh = ', ln_linssh 414 WRITE(numout,*) ' create mesh/mask file ln_meshmask = ', ln_meshmask 415 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' [m]' 416 WRITE(numout,*) ' ocean time step rn_dt = ', rn_dt 417 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 418 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs 419 ENDIF 420 ! 421 IF( TRIM(Agrif_CFixed()) == '0' ) THEN 435 ! 436 IF( TRIM(Agrif_CFixed()) == '0' ) THEN !set output file type for XIOS based on NEMO namelist 422 437 lrxios = ln_xios_read.AND.ln_rstart 423 !set output file type for XIOS based on NEMO namelist424 438 IF (nn_wxios > 0) lwxios = .TRUE. 425 439 nxioso = nn_wxios … … 553 567 CALL iom_getatt( inum, 'cn_cfg', cd_cfg ) ! returns ! if not found 554 568 CALL iom_getatt( inum, 'nn_cfg', kk_cfg ) ! returns -999 if not found 555 IF( TRIM(cd_cfg) == '!') cd_cfg = 'UNKNOWN'556 IF( kk_cfg == -999 ) kk_cfg = -9999999569 IF( TRIM(cd_cfg) == '!') cd_cfg = 'UNKNOWN' 570 IF( kk_cfg == -999 ) kk_cfg = -9999999 557 571 ENDIF 558 572 ! … … 701 715 ! 702 716 END SUBROUTINE cfg_write 703 717 704 718 !!====================================================================== 705 719 END MODULE domain -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dommsk.F90
r9657 r10001 202 202 ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 ) 203 203 ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 ) 204 ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 ) 204 205 205 206 -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90
r9939 r10001 4 4 !! Ocean : 5 5 !!====================================================================== 6 !! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code 7 !! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 6 !! History : 5.0 ! 2018-07 (G. Madec) Flux Form with Kinetic Energy conservation 7 !! ==>>> here z* and s* only (no z-tilde) 8 9 ! 1- remove z-tilde ==>>> pure z-star (or s-star) 10 ! 2- remove dom_vvl_interpol 11 10 12 !!---------------------------------------------------------------------- 11 13 … … 14 16 !! dom_vvl_sf_nxt : Compute next vertical scale factors 15 17 !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid 16 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another17 18 !! dom_vvl_rst : read/write restart file 18 19 !! dom_vvl_ctl : Check the vvl options … … 38 39 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 40 PUBLIC dom_vvl_sf_swp ! called by step.F90 40 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 PUBLIC ssh2e3_before ! ... 42 PUBLIC ssh2e3_now ! ... 41 43 42 44 ! !!* Namelist nam_vvl … … 61 63 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 62 64 65 66 !!gm add 67 !!gm 68 69 70 63 71 !! * Substitutions 64 72 # include "vectopt_loop_substitute.h90" … … 74 82 !! *** FUNCTION dom_vvl_alloc *** 75 83 !!---------------------------------------------------------------------- 84 ! 76 85 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 86 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN … … 115 124 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 125 !!---------------------------------------------------------------------- 117 INTEGER :: ji, jj, jk118 INTEGER :: ii0, ii1, ij0, ij1119 REAL(wp):: zcoef, z1_Dt120 126 !!---------------------------------------------------------------------- 121 127 ! … … 129 135 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 130 136 ! 131 ! ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf137 ! ! Read or initialize e3t_(b/n), ssh(b/n) 132 138 CALL dom_vvl_rst( nit000, 'READ' ) 133 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 134 ! 135 ! !== Set of all other vertical scale factors ==! (now and before) 136 ! ! Horizontal interpolation of e3t 137 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U 138 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 139 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V 140 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 141 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F 142 ! ! Vertical interpolation of e3t,u,v 143 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W 144 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' ) 145 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW 146 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 147 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW 148 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 149 150 ! 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 ! 155 ! !== 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_wp 158 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 159 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 160 gdepw_b(:,:,1) = 0.0_wp 161 DO jk = 2, jpk ! vertical sum 162 DO jj = 1,jpj 163 DO ji = 1,jpi 164 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 165 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 166 ! ! 0.5 where jk = mikt 167 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 168 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)) 176 END DO 177 END DO 178 END DO 179 ! 180 ! !== 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) 186 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) 192 END DO 193 ! 194 ! !== inverse of water column thickness ==! (u- and v- points) 195 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 196 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(:,:) ) 199 200 ! !== z_tilde coordinate case ==! (Restoring frequencies) 201 IF( ln_vvl_ztilde ) THEN 202 !!gm : idea: add here a READ in a file of custumized restoring frequency 203 ! ! Values in days provided via the namelist 204 ! ! use rsmall to avoid possible division by zero errors with faulty settings 205 frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 206 frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 207 ! 208 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 209 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 210 frq_rst_hdv(:,:) = 1._wp / rn_Dt 211 ENDIF 212 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 213 z1_Dt = 1._wp / rn_Dt 214 DO jj = 1, jpj 215 DO ji = 1, jpi 216 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 217 IF( ABS(gphit(ji,jj)) >= 6.) THEN 218 ! values outside the equatorial band and transition zone (ztilde) 219 frq_rst_e3t(ji,jj) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400._wp ) 220 frq_rst_hdv(ji,jj) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400._wp ) 221 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 222 ! values inside the equatorial band (ztilde as zstar) 223 frq_rst_e3t(ji,jj) = 0._wp 224 frq_rst_hdv(ji,jj) = z1_Dt 225 ELSE ! transition band (2.5 to 6 degrees N/S) 226 ! ! (linearly transition from z-tilde to z-star) 227 frq_rst_e3t(ji,jj) = 0._wp + ( frq_rst_e3t(ji,jj) - 0._wp ) * 0.5_wp & 228 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp ) ) 229 frq_rst_hdv(ji,jj) = z1_Dt + ( frq_rst_hdv(ji,jj) - z1_Dt ) * 0.5_wp & 230 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp ) ) 231 ENDIF 232 END DO 233 END DO 234 IF( cn_cfg == "orca" .AND. nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 235 ii0 = 103 ; ii1 = 111 236 ij0 = 128 ; ij1 = 135 ; 237 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp 238 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = z1_Dt 239 ENDIF 240 ENDIF 241 ENDIF 242 ! 243 IF(lwxios) THEN 244 ! define variables in restart file when writing with XIOS 139 ! 140 ! !== Set of all other vertical mesh fields ==! (now and before) 141 ! 142 ! !* BEFORE fields : 143 CALL ssh2e3_before ! set: hu , hv , r1_hu, r1_hv 144 ! ! e3t, e3w, e3u, e3uw, e3v, e3vw 145 ! 146 ! ! set one for all last level to the e3._0 value 147 e3t_b(:,:,jpk) = e3t_0(:,:,jpk) ; e3u_b(:,:,jpk) = e3w_0(:,:,jpk) ; e3v_b(:,:,jpk) = e3v_0(:,:,jpk) 148 e3w_b(:,:,jpk) = e3w_0(:,:,jpk) ; e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk) ; e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk) 149 ! 150 ! !* NOW fields : 151 CALL ssh2e3_now ! set: ht , hu , hv , r1_hu, r1_hv 152 ! ! e3t, e3w, e3u, e3uw, e3v, e3vw, e3f 153 ! ! gdept_n, gdepw_n, gde3w_n 154 ! 155 ! ! set one for all last level to the e3._0 value 156 e3t_n(:,:,jpk) = e3t_0(:,:,jpk) ; e3u_n(:,:,jpk) = e3w_0(:,:,jpk) ; e3v_n(:,:,jpk) = e3v_0(:,:,jpk) 157 e3w_n(:,:,jpk) = e3w_0(:,:,jpk) ; e3uw_n(:,:,jpk) = e3uw_0(:,:,jpk) ; e3vw_n(:,:,jpk) = e3vw_0(:,:,jpk) 158 e3f_n(:,:,jpk) = e3f_0(:,:,jpk) 159 ! 160 ! !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation) 161 e3t_a(:,:,:) = e3t_n(:,:,:) ; e3u_a(:,:,:) = e3u_n(:,:,:) ; e3v_a(:,:,:) = e3v_n(:,:,:) 162 163 !!gm 164 !!gm ===>>>> below: some issues to think about !!! 165 !!gm 166 !!gm fmask definition checked (0 or 1 nothing else) 167 !! in z-tilde or ALE e3._0 should be the time varying fields step forward with an euler scheme 168 !!gm e3w_b & gdept_b are not used except its update in agrif 169 !! and used to compute before slope of surface in dynldf_iso ==>>> remove it !!!! 170 !! NB: in triads on TRA, gdept_n is used !!!! BUG? 171 !!gm e3f_n almost not used ===>>>> verify whether it can be removed or not... 172 !! verify the use of wumask & wvmask mau be replaced by a multiplication by umask(k)*umask(k+1) ??? 173 !! 174 !!gm ISF case to be checked by Pierre Mathiot 175 !! 176 !!gm setting of e3._a for agrif.... TO BE CHECKED 177 178 ! 179 IF(lwxios) THEN ! define variables in restart file when writing with XIOS 245 180 CALL iom_set_rstw_var_active('e3t_b') 246 181 CALL iom_set_rstw_var_active('e3t_n') 247 ! ! ----------------------- !248 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !249 ! ! ----------------------- !250 CALL iom_set_rstw_var_active('tilde_e3t_b')251 CALL iom_set_rstw_var_active('tilde_e3t_n')252 END IF253 ! ! -------------!254 IF( ln_vvl_ztilde ) THEN ! z_tilde case !255 ! ! ------------ !256 CALL iom_set_rstw_var_active('hdiv_lf')257 ENDIF258 !259 182 ENDIF 260 183 ! … … 287 210 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 288 211 ! 289 INTEGER :: ji, jj, jk ! dummy loop indices 290 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 291 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars 292 LOGICAL :: ll_do_bclinic ! local logical 293 REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv 294 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t 212 INTEGER :: ji, jj, jk ! dummy loop indices 213 REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h 295 214 !!---------------------------------------------------------------------- 296 215 ! … … 305 224 ENDIF 306 225 307 ll_do_bclinic = .TRUE. 308 IF( PRESENT(kcall) ) THEN 309 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 310 ENDIF 311 312 ! ******************************* ! 313 ! After acale factors at t-points ! 314 ! ******************************* ! 315 ! ! --------------------------------------------- ! 316 ! ! z_star coordinate and barotropic z-tilde part ! 317 ! ! --------------------------------------------- ! 318 ! 319 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 320 DO jk = 1, jpkm1 321 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 322 e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 323 END DO 324 ! 325 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 326 ! ! ------baroclinic part------ ! 327 ! I - initialization 328 ! ================== 329 330 ! 1 - barotropic divergence 331 ! ------------------------- 332 zhdiv(:,:) = 0._wp 333 zht(:,:) = 0._wp 334 DO jk = 1, jpkm1 335 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 336 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 337 END DO 338 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 339 340 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) 341 ! -------------------------------------------------- 342 IF( ln_vvl_ztilde ) THEN 343 IF( kt > nit000 ) THEN 344 DO jk = 1, jpkm1 345 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:) & 346 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 347 END DO 348 ENDIF 349 ENDIF 350 351 ! II - after z_tilde increments of vertical scale factors 352 ! ======================================================= 353 te3t_a(:,:,:) = 0._wp ! te3t_a used to store tendency terms 354 355 ! 1 - High frequency divergence term 356 ! ---------------------------------- 357 IF( ln_vvl_ztilde ) THEN ! z_tilde case 358 DO jk = 1, jpkm1 359 te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 360 END DO 361 ELSE ! layer case 362 DO jk = 1, jpkm1 363 te3t_a(:,:,jk) = te3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 364 END DO 365 ENDIF 366 367 ! 2 - Restoring term (z-tilde case only) 368 ! ------------------ 369 IF( ln_vvl_ztilde ) THEN 370 DO jk = 1, jpk 371 te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk) 372 END DO 373 ENDIF 374 375 ! 3 - Thickness diffusion term 376 ! ---------------------------- 377 zwu(:,:) = 0._wp 378 zwv(:,:) = 0._wp 379 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 380 DO jj = 1, jpjm1 381 DO ji = 1, fs_jpim1 ! vector opt. 382 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 383 & * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj ,jk) ) 384 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 385 & * ( te3t_b(ji,jj,jk) - te3t_b(ji ,jj+1,jk) ) 386 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 387 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 388 END DO 389 END DO 390 END DO 391 DO jj = 1, jpj ! b - correction for last oceanic u-v points 392 DO ji = 1, jpi 393 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 394 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 395 END DO 396 END DO 397 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 398 DO jj = 2, jpjm1 399 DO ji = fs_2, fs_jpim1 ! vector opt. 400 te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 401 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 402 & ) * r1_e1e2t(ji,jj) 403 END DO 404 END DO 405 END DO 406 ! ! d - thickness diffusion transport: boundary conditions 407 ! (stored for tracer advction and continuity equation) 408 CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp) 409 410 ! 4 - Time stepping of baroclinic scale factors 411 ! --------------------------------------------- 412 ! Leapfrog time stepping 413 ! ~~~~~~~~~~~~~~~~~~~~~~ 414 CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp ) 415 te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:) 416 417 ! Maximum deformation control 418 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 419 ze3t(:,:,jpk) = 0._wp 420 DO jk = 1, jpkm1 421 ze3t(:,:,jk) = te3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 422 END DO 423 z_tmax = MAXVAL( ze3t(:,:,:) ) 424 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 425 z_tmin = MINVAL( ze3t(:,:,:) ) 426 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 427 ! - ML - test: for the moment, stop simulation for too large e3_t variations 428 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 429 IF( lk_mpp ) THEN 430 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 431 CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 432 ELSE 433 ijk_max = MAXLOC( ze3t(:,:,:) ) 434 ijk_max(1) = ijk_max(1) + nimpp - 1 435 ijk_max(2) = ijk_max(2) + njmpp - 1 436 ijk_min = MINLOC( ze3t(:,:,:) ) 437 ijk_min(1) = ijk_min(1) + nimpp - 1 438 ijk_min(2) = ijk_min(2) + njmpp - 1 439 ENDIF 440 IF (lwp) THEN 441 WRITE(numout, *) 'MAX( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 442 WRITE(numout, *) 'at i, j, k=', ijk_max 443 WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 444 WRITE(numout, *) 'at i, j, k=', ijk_min 445 CALL ctl_warn('MAX( ABS( te3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 446 ENDIF 447 ENDIF 448 ! - ML - end test 449 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 450 te3t_a(:,:,:) = MIN( te3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) ) 451 te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 452 453 ! 454 ! "tilda" change in the after scale factor 455 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 456 DO jk = 1, jpkm1 457 dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_b(:,:,jk) 458 END DO 459 ! III - Barotropic repartition of the sea surface height over the baroclinic profile 460 ! ================================================================================== 461 ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n) 462 ! - ML - baroclinicity error should be better treated in the future 463 ! i.e. locally and not spread over the water column. 464 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 465 zht(:,:) = 0._wp 466 DO jk = 1, jpkm1 467 zht(:,:) = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk) 468 END DO 469 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 470 DO jk = 1, jpkm1 471 dte3t_a(:,:,jk) = dte3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 472 END DO 473 474 ENDIF 475 476 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 477 ! ! ---baroclinic part--------- ! 478 DO jk = 1, jpkm1 479 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dte3t_a(:,:,jk) * tmask(:,:,jk) 480 END DO 481 ENDIF 482 483 IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging 484 ! 485 IF( lwp ) WRITE(numout, *) 'kt =', kt 486 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 487 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 488 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 489 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(te3t_a))) =', z_tmax 490 END IF 491 ! 492 zht(:,:) = 0.0_wp 493 DO jk = 1, jpkm1 494 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 495 END DO 496 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 497 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 498 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 499 ! 500 zht(:,:) = 0.0_wp 501 DO jk = 1, jpkm1 502 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 503 END DO 504 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 505 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 506 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 507 ! 508 zht(:,:) = 0.0_wp 509 DO jk = 1, jpkm1 510 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 511 END DO 512 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 513 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 514 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 515 ! 516 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) 517 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 518 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 519 ! 520 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshn(:,:) ) ) 521 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 522 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 523 ! 524 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssha(:,:) ) ) 525 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 526 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 527 END IF 528 529 ! *********************************** ! 530 ! After scale factors at u- v- points ! 531 ! *********************************** ! 532 533 CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 534 CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 535 536 ! *********************************** ! 537 ! After depths at u- v points ! 538 ! *********************************** ! 539 540 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 541 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 542 DO jk = 2, jpkm1 543 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 544 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 545 END DO 546 ! ! Inverse of the local depth 547 !!gm BUG ? don't understand the use of umask_i here ..... 226 ! ! ------------------! 227 ! ! z_star coordinate ! (and barotropic z-tilde part) 228 ! ! ------------------! 229 ! 230 ! !== after ssh ==! (u- and v-points) 231 DO jj = 2, jpjm1 ; DO ji = 2, jpim1 232 zsshu_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji+1,jj) ) * ssumask(ji,jj) 233 zsshv_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji,jj+1) ) * ssvmask(ji,jj) 234 END DO ; END DO 235 CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp ) 236 ! 237 ! !== after depths and its inverse ==! 238 hu_a(:,:) = hu_0(:,:) + zsshu_h(:,:) 239 hv_a(:,:) = hv_0(:,:) + zsshv_h(:,:) 548 240 r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 549 241 r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 242 ! 243 ! !== after scale factors ==! (e3t , e3u , e3v) 244 zssht_h(:,:) = ssha (:,:) * r1_ht_0(:,:) ! t-point 245 zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! u-point 246 zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! v-point 247 DO jk = 1, jpkm1 248 e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 249 e3u_a(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) 250 e3v_a(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) 251 END DO 550 252 ! 551 253 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 581 283 ! 582 284 INTEGER :: ji, jj, jk ! dummy loop indices 583 REAL(wp) :: zcoef, ze3f ! local scalar285 REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h, zsshf_h 584 286 !!---------------------------------------------------------------------- 585 287 ! … … 594 296 ENDIF 595 297 ! 596 ! Time filter and swap of scale factors 597 ! ===================================== 598 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 599 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 600 IF( l_1st_euler ) THEN 601 te3t_n(:,:,:) = te3t_a(:,:,:) 602 ELSE 603 DO jk = 1, jpk 604 DO jj = 1, jpj 605 DO ji = 1, jpi 606 ze3f = te3t_n(ji,jj,jk) & 607 & + rn_atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) ) 608 te3t_b(ji,jj,jk) = ze3f 609 te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk) 610 END DO 611 END DO 612 END DO 613 ENDIF 614 ENDIF 615 gdept_b(:,:,:) = gdept_n(:,:,:) 616 gdepw_b(:,:,:) = gdepw_n(:,:,:) 617 618 e3t_n(:,:,:) = e3t_a(:,:,:) 619 e3u_n(:,:,:) = e3u_a(:,:,:) 620 e3v_n(:,:,:) = e3v_a(:,:,:) 621 622 ! Compute all missing vertical scale factor and depths 623 ! ==================================================== 624 ! Horizontal scale factor interpolations 625 ! -------------------------------------- 298 ! Swap and Compute all missing vertical scale factor and depths 299 ! ============================================================= 626 300 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 627 301 ! - JC - hu_b, hv_b, hur_b, hvr_b also 628 629 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 630 631 ! Vertical scale factor interpolations 632 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 633 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 634 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 635 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' ) 636 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 637 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 638 639 ! t- and w- points depth (set the isf depth as it is in the initial step) 640 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 641 gdepw_n(:,:,1) = 0.0_wp 642 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 643 DO jk = 2, jpk 644 DO jj = 1,jpj 645 DO ji = 1,jpi 646 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 647 ! 1 for jk = mikt 648 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 649 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 650 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 651 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 652 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 653 END DO 654 END DO 655 END DO 656 302 ! 303 ! - GM - to be updated : e3f_n, e3w_n , e3uw_n , e3vw_n 304 ! e3w_b , e3uw_b , e3vw_b 305 ! gdept_n , gdepw_n , gde3w_n 306 ! ht_n 307 ! to be swap : hu_a , hv_a , r1_hu_a , r1_hv_a 308 ! 657 309 ! Local depth and Inverse of the local depth of the water 658 310 ! ------------------------------------------------------- 311 ! ! swap of depth and scale factors 312 ! ! =============================== 313 DO jk = 1, jpkm1 ! swap n--> a 314 gdept_b(:,:,jk) = gdept_n(:,:,jk) ! depth at t and w 315 gdepw_b(:,:,jk) = gdepw_n(:,:,jk) 316 e3t_n (:,:,jk) = e3t_a (:,:,jk) ! e3t, e3u, e3v 317 e3u_n (:,:,jk) = e3u_a (:,:,jk) 318 e3v_n (:,:,jk) = e3v_a (:,:,jk) 319 END DO 320 ht_n(:,:) = ht_0(:,:) + sshn(:,:) ! ocean thickness 321 ! 659 322 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 660 323 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 661 324 ! 662 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 663 DO jk = 2, jpkm1 664 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 325 ! !== before : 326 ! !* ssh at u- and v-points) 327 DO jj = 2, jpjm1 ; DO ji = 2, jpim1 328 zsshu_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) 329 zsshv_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) 330 END DO ; END DO 331 CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp ) 332 ! 333 ! !* e3w_b , e3uw_b , e3vw_b 334 zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:) ! w-point 335 zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! uw-point 336 zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! vw-point 337 DO jk = 1, jpkm1 338 e3w_b(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 339 e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) 340 e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) 665 341 END DO 666 342 ! 343 ! !== now : 344 ! !* ssh at u- and v-points) 345 DO jj = 1, jpjm1 ; DO ji = 1, jpim1 ! start from 1 for f-point 346 zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji+1,jj ) ) * ssumask(ji,jj) 347 zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) ) * ssvmask(ji,jj) 348 zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 349 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 350 END DO ; END DO 351 CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp ) 352 ! 353 ! !* e3w_n , e3uw_n , e3vw_n, e3f_n 354 zssht_h(:,:) = sshn (:,:) * r1_ht_0(:,:) ! t- & w-point 355 zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! uw-point 356 zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! vw-point 357 zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:) ! f-point 358 DO jk = 1, jpkm1 359 e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 360 e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) ) 361 e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) ) 362 e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) ) 363 END DO 364 ! 365 zssht_h(:,:) = 1._wp + sshn (:,:) * r1_ht_0(:,:) ! t-point 366 ! 367 IF( ln_isfcav ) THEN ! ISF cavities : ssh scaling not applied over the iceshelf thickness 368 DO jk = 1, jpkm1 369 gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 370 gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 371 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn (:,:) 372 END DO 373 ELSE ! no ISF cavities 374 DO jk = 1, jpkm1 375 gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:) 376 gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:) 377 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 378 END DO 379 ENDIF 380 ! 667 381 ! write restart file 668 382 ! ================== … … 672 386 ! 673 387 END SUBROUTINE dom_vvl_sf_swp 674 675 676 SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )677 !!---------------------------------------------------------------------678 !! *** ROUTINE dom_vvl__interpol ***679 !!680 !! ** Purpose : interpolate scale factors from one grid point to another681 !!682 !! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)683 !! - horizontal interpolation: grid cell surface averaging684 !! - vertical interpolation: simple averaging685 !!----------------------------------------------------------------------686 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated687 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3688 CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors689 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'690 !691 INTEGER :: ji, jj, jk ! dummy loop indices692 REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F693 !!----------------------------------------------------------------------694 !695 IF(ln_wd_il) THEN696 zlnwd = 1.0_wp697 ELSE698 zlnwd = 0.0_wp699 END IF700 !701 SELECT CASE ( pout ) !== type of interpolation ==!702 !703 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean704 DO jk = 1, jpk705 DO jj = 1, jpjm1706 DO ji = 1, fs_jpim1 ! vector opt.707 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &708 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &709 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )710 END DO711 END DO712 END DO713 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )714 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)715 !716 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean717 DO jk = 1, jpk718 DO jj = 1, jpjm1719 DO ji = 1, fs_jpim1 ! vector opt.720 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &721 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &722 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )723 END DO724 END DO725 END DO726 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )727 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)728 !729 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean730 DO jk = 1, jpk731 DO jj = 1, jpjm1732 DO ji = 1, fs_jpim1 ! vector opt.733 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &734 & * r1_e1e2f(ji,jj) &735 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &736 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )737 END DO738 END DO739 END DO740 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )741 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)742 !743 CASE( 'W' ) !* from T- to W-point : vertical simple mean744 !745 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)746 ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing747 !!gm BUG? use here wmask in case of ISF ? to be checked748 DO jk = 2, jpk749 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &750 & * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &751 & + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &752 & * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )753 END DO754 !755 CASE( 'UW' ) !* from U- to UW-point : vertical simple mean756 !757 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)758 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing759 !!gm BUG? use here wumask in case of ISF ? to be checked760 DO jk = 2, jpk761 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &762 & * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &763 & + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &764 & * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )765 END DO766 !767 CASE( 'VW' ) !* from V- to VW-point : vertical simple mean768 !769 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)770 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing771 !!gm BUG? use here wvmask in case of ISF ? to be checked772 DO jk = 2, jpk773 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &774 & * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &775 & + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &776 & * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )777 END DO778 END SELECT779 !780 END SUBROUTINE dom_vvl_interpol781 388 782 389 … … 804 411 IF( ln_rstart ) THEN !* Read the restart file 805 412 CALL rst_read_open ! open the restart file if necessary 806 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios )807 413 ! 808 id1 = iom_varid( numror, 'e3t_b' , ldstop = .FALSE. ) 809 id2 = iom_varid( numror, 'e3t_n' , ldstop = .FALSE. ) 810 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 811 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 812 id5 = iom_varid( numror, 'hdiv_lf' , ldstop = .FALSE. ) 813 ! ! --------- ! 814 ! ! all cases ! 815 ! ! --------- ! 414 id1 = iom_varid( numror, 'sshb' , ldstop = .FALSE. ) 415 id2 = iom_varid( numror, 'sshn' , ldstop = .FALSE. ) 416 ! 816 417 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 817 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 818 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 819 ! needed to restart if land processor not computed 820 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 821 WHERE ( tmask(:,:,:) == 0.0_wp ) 822 e3t_n(:,:,:) = e3t_0(:,:,:) 823 e3t_b(:,:,:) = e3t_0(:,:,:) 824 END WHERE 418 IF(lwp) write(numout,*) 'dom_vvl_rst : both sshb and sshn found in restart files' 419 ! 420 !!gm Question: use jpdom_data above to read data over jpi x jpj (like is dom_hgr_read and dom_zgr_read) 421 !! so that it will work with land processor suppression 422 ! CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios ) 423 ! CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb, ldxios = lrxios ) 424 !!gm 425 CALL iom_get( numror, jpdom_data, 'sshn' , sshn, ldxios = lrxios ) 426 CALL iom_get( numror, jpdom_data, 'sshb' , sshb, ldxios = lrxios ) 427 !!gm end 825 428 IF( l_1st_euler ) THEN 826 e3t_b(:,:,:) = e3t_n(:,:,:)429 sshb(:,:) = sshn(:,:) 827 430 ENDIF 828 431 ELSE IF( id1 > 0 ) THEN 829 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'830 IF(lwp) write(numout,*) ' e3t_n set equal to e3t_b.'831 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 832 CALL iom_get( numror, jpdom_ autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )833 e3t_n(:,:,:) = e3t_b(:,:,:)432 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files' 433 IF(lwp) write(numout,*) ' set sshn = sshb and force l_1st_euler = true' 434 !!gm CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 435 CALL iom_get( numror, jpdom_data, 'sshb', sshb, ldxios = lrxios ) 436 sshn(:,:) = sshb(:,:) 834 437 l_1st_euler = .TRUE. 835 438 ELSE IF( id2 > 0 ) THEN 836 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 837 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 838 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 839 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 840 e3t_b(:,:,:) = e3t_n(:,:,:) 439 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files' 440 IF(lwp) write(numout,*) 'set sshb = sshn and force l_1st_euler = true' 441 CALL iom_get( numror, jpdom_data, 'sshn', sshb, ldxios = lrxios ) 442 sshb(:,:) = sshn(:,:) 841 443 l_1st_euler = .TRUE. 842 444 ELSE 843 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 844 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 845 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 846 DO jk = 1, jpk 847 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 848 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 849 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 850 END DO 851 e3t_b(:,:,:) = e3t_n(:,:,:) 445 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file' 446 IF(lwp) write(numout,*) 'set sshb = sshn = 0 and force l_1st_euler = true' 447 sshb(:,:) = 0._wp 448 sshn(:,:) = 0._wp 852 449 l_1st_euler = .TRUE. 853 450 ENDIF 854 ! ! ----------- !855 IF( ln_vvl_zstar ) THEN ! z_star case !856 ! ! ----------- !857 IF( MIN( id3, id4 ) > 0 ) THEN858 CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )859 ENDIF860 ! ! ----------------------- !861 ELSE ! z_tilde and layer cases !862 ! ! ----------------------- !863 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist864 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lrxios )865 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios )866 ELSE ! one at least array is missing867 te3t_b(:,:,:) = 0.0_wp868 te3t_n(:,:,:) = 0.0_wp869 ENDIF870 ! ! ------------ !871 IF( ln_vvl_ztilde ) THEN ! z_tilde case !872 ! ! ------------ !873 IF( id5 > 0 ) THEN ! required array exists874 CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )875 ELSE ! array is missing876 hdiv_lf(:,:,:) = 0.0_wp877 ENDIF878 ENDIF879 ENDIF880 !881 451 ELSE !* Initialize at "rest" 882 452 ! 883 884 453 IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential 885 454 ! 886 IF( cn_cfg == 'wad' ) THEN 887 ! Wetting and drying test case 455 IF( cn_cfg == 'wad' ) THEN ! Wetting and drying test case 888 456 CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb ) 889 tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones457 tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones 890 458 sshn (:,:) = sshb(:,:) 891 459 un (:,:,:) = ub (:,:,:) 892 460 vn (:,:,:) = vb (:,:,:) 893 ELSE 894 ! if not test case 461 ELSE ! Not the test case 895 462 sshn(:,:) = -ssh_ref 896 463 sshb(:,:) = -ssh_ref 897 464 ! 898 465 DO jj = 1, jpj 899 466 DO ji = 1, jpi 900 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 901 467 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 902 468 sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 903 469 sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 904 470 ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 905 471 ENDIF 906 ENDDO 907 ENDDO 908 ENDIF !If test case else 909 910 ! Adjust vertical metrics for all wad 911 DO jk = 1, jpk 912 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 913 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 914 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 915 END DO 916 e3t_b(:,:,:) = e3t_n(:,:,:) 917 472 END DO 473 END DO 474 ENDIF ! If test case else 475 ! 918 476 DO ji = 1, jpi 919 477 DO jj = 1, jpj 920 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN478 IF ( ht_0(ji,jj) /= 0._wp .AND. NINT( ssmask(ji,jj) ) == 1 ) THEN 921 479 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 922 480 ENDIF … … 926 484 ELSE 927 485 ! 928 ! Just to read set ssh in fact, called latter once vertical grid 929 ! is set up: 486 ! Just to read set ssh in fact, called latter once vertical grid is set up: 930 487 ! CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb ) 931 ! ! 932 ! DO jk=1,jpk 933 ! e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) & 934 ! & / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 935 ! END DO 936 ! e3t_n(:,:,:) = e3t_b(:,:,:) 937 sshn(:,:)=0._wp 938 e3t_n(:,:,:)=e3t_0(:,:,:) 939 e3t_b(:,:,:)=e3t_0(:,:,:) 488 sshn(:,:) = 0._wp 489 sshb(:,:) = 0._wp 940 490 ! 941 END IF ! end of ll_wd edits942 943 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN944 te3t_b(:,:,:) = 0._wp945 te3t_n(:,:,:) = 0._wp946 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp947 491 END IF 492 ! 948 493 ENDIF 949 494 ! 950 495 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 951 496 ! ! =================== 952 IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 953 IF( lwxios ) CALL iom_swap( cwxios_context ) 954 ! ! --------- ! 955 ! ! all cases ! 956 ! ! --------- ! 957 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios ) 958 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios ) 959 ! ! ----------------------- ! 960 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! 961 ! ! ----------------------- ! 962 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios) 963 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios) 964 END IF 965 ! ! -------------! 966 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 967 ! ! ------------ ! 968 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios) 969 ENDIF 970 ! 971 IF( lwxios ) CALL iom_swap( cxios_context ) 497 498 !!gm DO NOTHING, sshb and sshn are written in restart.F90 499 972 500 ENDIF 973 501 ! … … 1024 552 ENDIF 1025 553 ! 554 555 !!gm 556 IF ( ln_vvl_ztilde .OR. ln_vvl_ztilde_as_zstar ) CALL ctl_stop( 'z-tilde NOT available in this branch' ) 557 !!gm 558 559 ! 1026 560 ioptio = 0 ! Parameter control 1027 561 IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true. … … 1047 581 END SUBROUTINE dom_vvl_ctl 1048 582 583 584 SUBROUTINE ssh2e3_now 585 !!---------------------------------------------------------------------- 586 !! *** ROUTINE ssh2e3_now *** 587 !!---------------------------------------------------------------------- 588 INTEGER :: ji, jj, jk 589 REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h, zsshf_h 590 !!---------------------------------------------------------------------- 591 ! 592 ! !== ssh at u- and v-points ==! 593 ! 594 DO jj = 1, jpjm1 ! start from 1 due to f-point 595 DO ji = 1, jpim1 596 zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji+1,jj ) ) * ssumask(ji,jj) 597 zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) ) * ssvmask(ji,jj) 598 zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 599 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 600 END DO 601 END DO 602 CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp ) 603 ! 604 ! !== ht, hu and hv == ! (and their inverse) 605 ! 606 ht_n (:,:) = ht_0(:,:) + sshn (:,:) 607 hu_n (:,:) = hu_0(:,:) + zsshu_h(:,:) 608 hv_n (:,:) = hv_0(:,:) + zsshv_h(:,:) 609 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) ! ss mask mask due to ISF 610 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 611 ! 612 ! !== ssh / h factor at t-, u- ,v- & f-points ==! 613 ! 614 zssht_h(:,:) = sshn (:,:) * r1_ht_0(:,:) 615 zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) 616 zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) 617 zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:) 618 ! 619 ! !== e3t, e3w , e3u, e3uw , e3v, e3vw , and e3f ==! 620 ! 621 DO jk = 1, jpkm1 622 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 623 e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 624 ! 625 e3u_n(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) 626 e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) ) 627 ! 628 e3v_n(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) 629 e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) ) 630 ! 631 e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) ) 632 END DO 633 ! 634 ! !== depth of t- and w-points ==! 635 ! 636 zssht_h(:,:) = 1._wp + zssht_h(:,:) ! = 1 + sshn / ht_0 637 ! 638 IF( ln_isfcav ) THEN ! ISF cavities : ssh scaling not applied over the iceshelf thickness 639 DO jk = 1, jpkm1 640 gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 641 gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 642 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 643 END DO 644 ELSE ! no ISF cavities 645 DO jk = 1, jpkm1 646 gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:) 647 gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:) 648 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 649 END DO 650 ENDIF 651 ! 652 END SUBROUTINE ssh2e3_now 653 654 655 SUBROUTINE ssh2e3_before 656 !!---------------------------------------------------------------------- 657 !! *** ROUTINE ssh2e3_before *** 658 !!---------------------------------------------------------------------- 659 INTEGER :: ji, jj, jk 660 REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h 661 !!---------------------------------------------------------------------- 662 ! 663 ! !== ssh at u- and v-points ==! 664 DO jj = 2, jpjm1 665 DO ji = 2, jpim1 666 zsshu_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) 667 zsshv_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) 668 END DO 669 END DO 670 CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp ) 671 ! 672 ! !== ht, hu and hv == ! (and their inverse) 673 hu_b(:,:) = hu_0(:,:) + zsshu_h(:,:) 674 hv_b(:,:) = hv_0(:,:) + zsshv_h(:,:) 675 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 676 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 677 ! 678 ! 679 ! !== ssh / h factor at t-, u- ,v- & f-points ==! 680 zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:) 681 zsshu_h (:,:) = zsshu_h(:,:) * r1_hu_0(:,:) 682 zsshv_h (:,:) = zsshv_h(:,:) * r1_hv_0(:,:) 683 ! 684 ! !== e3t, e3w , e3u, e3uw , and e3v, e3vw ==! 685 DO jk = 1, jpkm1 686 e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 687 e3w_b(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 688 ! 689 e3u_b(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h (:,:) * umask(:,:,jk) ) 690 e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h (:,:) * wumask(:,:,jk) ) 691 ! 692 e3v_b(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h (:,:) * vmask(:,:,jk) ) 693 e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h (:,:) * wvmask(:,:,jk) ) 694 END DO 695 ! 696 END SUBROUTINE ssh2e3_before 697 1049 698 !!====================================================================== 1050 699 END MODULE domvvl -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90
r9939 r10001 14 14 USE dom_oce ! ocean space and time domain 15 15 USE domwri ! ocean space and time domain 16 USE domvvl , ONLY : dom_vvl_interpol17 16 USE phycst ! physical constants 18 17 USE sbc_oce ! surface boundary condition variables … … 137 136 REAL(wp), DIMENSION(jpi,jpj) :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t 138 137 REAL(wp), DIMENSION(jpi,jpj) :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1 139 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn , ztrp138 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn 140 139 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask1, zwmaskb, ztmp3d 141 140 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0 141 REAL(wp), DIMENSION(jpi,jpj) :: z_ssh_h0, zsshu, zsshv, zsshf 142 142 !!---------------------------------------------------------------------- 143 143 ! … … 163 163 DO jj = 2,jpj-1 164 164 DO ji = fs_2, fs_jpim1 ! vector opt. 165 jip1=ji+1; jim1=ji-1; 166 jjp1=jj+1; jjm1=jj-1; 167 summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1)) 165 summsk = zsmask0(ji+1,jj)+zsmask0(ji-1,jj)+zsmask0(ji,jj+1)+zsmask0(ji,jj-1) 168 166 IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 169 sshn(ji,jj)=( zssh0(ji p1,jj)*zsmask0(jip1,jj) &170 & + zssh0(ji m1,jj)*zsmask0(jim1,jj) &171 & + zssh0(ji,jj p1)*zsmask0(ji,jjp1) &172 & + zssh0(ji,jj m1)*zsmask0(ji,jjm1))/summsk173 zsmask1(ji,jj) =1._wp167 sshn(ji,jj)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj) & 168 & + zssh0(ji-1,jj)*zsmask0(ji-1,jj) & 169 & + zssh0(ji,jj+1)*zsmask0(ji,jj+1) & 170 & + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk 171 zsmask1(ji,jj) = 1._wp 174 172 ENDIF 175 173 END DO … … 185 183 IF ( .NOT.ln_linssh ) THEN 186 184 ! Reconstruction of all vertical scale factors at now time steps 187 ! ============================================================================= 185 ! ====================================================================== 186 187 !!gm Question : bug ???? 188 ! at this stage is ht_0, r1_ht0 already computed ???? 189 ! I have the feeling that this should be done in a different manner... 190 ! that is iscpl_rst_interpol should be call directly in restart !!! 191 192 ! in my changes here I assume that ht_0 AND r1_ht_0 have been updated 193 ! Note that the former calculation were using ht_0 so if it as not been updated ===>>> BUG 194 !!gm 195 196 188 197 ! Horizontal scale factor interpolations 189 198 ! -------------------------------------- 190 DO jk = 1,jpk 191 DO jj=1,jpj 192 DO ji=1,jpi 193 IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN 194 e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / & 195 & ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) ) 196 ENDIF 197 END DO 198 END DO 199 END DO 200 ! 201 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 202 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 203 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 204 205 ! Vertical scale factor interpolations 206 ! ------------------------------------ 207 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 208 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 209 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 IF ( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp ) THEN 202 DO jk = 1, jpk 203 e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) ) 204 END DO 205 ENDIF 206 END DO 207 END DO 208 ! 209 !!gm Note that if this routine is called in dom_vvl_init then all the lines below are uselss !!! 210 !! they are a duplication of dom_vvl_init lines 211 212 ! !== now fields ==! 213 ! 214 ! !* ssh at u- and v-points) 215 DO jj = 1, jpjm1 ; DO ji = 1, jpim1 ! start from 1 due to f-point 216 zsshu(ji,jj) = 0.5_wp * ( sshn(ji ,jj) + sshn(ji+1,jj ) ) * ssumask(ji,jj) 217 zsshv(ji,jj) = 0.5_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) ) * ssvmask(ji,jj) 218 zsshf(ji,jj) = 0.25_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 219 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 220 END DO ; END DO 221 CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp , zsshf(:,:),'F', 1._wp ) 222 ! 223 ! !* hu and hv (and their inverse) 224 ht_n (:,:) = ht_0(:,:) + sshn(:,:) 225 hu_n (:,:) = hu_0(:,:) + zsshu(:,:) 226 hv_n (:,:) = hv_0(:,:) + zsshv(:,:) 227 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) ! ss mask mask due to ISF 228 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 229 ! 230 ! !* e3u, e3uw and e3v, e3vw 231 z_ssh_h0(:,:) = sshn (:,:) * r1_ht_0(:,:) ! t-point 232 DO jk = 1, jpkm1 233 e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * tmask(:,:,jk) ) 234 END DO 235 z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:) ! u-point 236 DO jk = 1, jpkm1 237 e3u_n (:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * umask(:,:,jk) ) 238 e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wumask(:,:,jk) ) 239 END DO 240 z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:) ! v-point 241 DO jk = 1, jpkm1 242 e3v_n (:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) ) 243 e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wvmask(:,:,jk) ) 244 END DO 245 z_ssh_h0(:,:) = zsshf(:,:) * r1_hf_0(:,:) ! f-point 246 DO jk = 1, jpkm1 247 e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * fmask(:,:,jk) ) 248 END DO 249 250 z_ssh_h0(:,:) = 1._wp + sshn(:,:) * r1_ht_0(:,:) ! t-point 251 ! 252 IF( ln_isfcav ) THEN ! iceshelf cavities : ssh scaling not applied over the iceshelf thickness 253 DO jk = 1, jpkm1 254 gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:) 255 gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:) 256 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 257 END DO 258 ELSE 259 DO jk = 1, jpkm1 260 gdept_n(:,:,jk) = gdept_0(:,:,jk) * z_ssh_h0(:,:) 261 gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * z_ssh_h0(:,:) 262 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 263 END DO 264 ENDIF 210 265 211 ! t- and w- points depth 212 ! ---------------------- 213 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 214 gdepw_n(:,:,1) = 0.0_wp 215 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 216 DO jj = 1,jpj 217 DO ji = 1,jpi 218 DO jk = 2,mikt(ji,jj)-1 219 gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 220 gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 221 gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 222 END DO 223 IF (mikt(ji,jj) > 1) THEN 224 jk = mikt(ji,jj) 225 gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk) 226 gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 227 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk ) - sshn (ji,jj) 228 END IF 229 DO jk = mikt(ji,jj)+1, jpk 230 gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) 231 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 232 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk ) - sshn (ji,jj) 233 END DO 234 END DO 235 END DO 236 237 ! t-, u- and v- water column thickness 238 ! ------------------------------------ 239 ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp 240 DO jk = 1, jpkm1 241 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 242 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 243 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 244 END DO 245 ! ! Inverse of the local depth 246 r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 247 r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 248 249 END IF 266 ENDIF 250 267 251 268 !============================================================================= 252 269 ! compute velocity 253 270 ! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor). 254 ub(:,:,:) =un(:,:,:)255 vb(:,:,:) =vn(:,:,:)256 DO jk = 1, jpk257 DO jj = 1, jpj258 DO ji = 1, jpi271 ub(:,:,:) = un(:,:,:) 272 vb(:,:,:) = vn(:,:,:) 273 DO jk = 1, jpk 274 DO jj = 1, jpj 275 DO ji = 1, jpi 259 276 un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk); 260 277 vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk); … … 265 282 ! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column) 266 283 ! compute barotropic velocity now and after 267 ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 268 zbub(:,:) = SUM(ztrp,DIM=3) 269 ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 270 zbvb(:,:) = SUM(ztrp,DIM=3) 271 ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 272 zbun(:,:) = SUM(ztrp,DIM=3) 273 ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 274 zbvn(:,:) = SUM(ztrp,DIM=3) 284 zbub(:,:) = SUM( ub(:,:,:)*pe3u_b(:,:,:) , DIM=3 ) 285 zbvb(:,:) = SUM( vb(:,:,:)*pe3v_b(:,:,:) , DIM=3 ) 286 zbun(:,:) = SUM( un(:,:,:)* e3u_n(:,:,:) , DIM=3 ) 287 zbvn(:,:) = SUM( vn(:,:,:)* e3v_n(:,:,:) , DIM=3 ) 275 288 276 289 ! new water column … … 285 298 zucorr = 0._wp 286 299 zvcorr = 0._wp 287 DO jj = 1, jpj288 DO ji = 1, jpi300 DO jj = 1, jpj 301 DO ji = 1, jpi 289 302 IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN 290 303 zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj) … … 297 310 298 311 ! update velocity 299 DO jk = 1, jpk300 un(:,:,jk) =(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)301 vn(:,:,jk) =(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)312 DO jk = 1, jpk 313 un(:,:,jk) = ( un(:,:,jk) - zucorr(:,:) ) * umask(:,:,jk) 314 vn(:,:,jk) = ( vn(:,:,jk) - zvcorr(:,:) ) * vmask(:,:,jk) 302 315 END DO 303 316 … … 305 318 ! compute temp and salt 306 319 ! compute new tn and sn if we open a new cell 307 tsb (:,:,:,:) = tsn(:,:,:,:)308 zts0 (:,:,:,:) = tsn(:,:,:,:)320 tsb (:,:,:,:) = tsn (:,:,:,:) 321 zts0 (:,:,:,:) = tsn (:,:,:,:) 309 322 ztmask1(:,:,:) = ptmask_b(:,:,:) 310 323 ztmask0(:,:,:) = ptmask_b(:,:,:) 311 324 DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case) 312 313 314 315 316 317 318 319 320 321 322 & 323 & 324 & 325 326 & 327 & 328 & 329 330 331 332 333 334 335 336 & +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk337 338 & +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk339 340 ENDIF341 ENDIF342 343 344 345 346 347 348 349 350 351 325 DO jk = 1,jpk-1 326 zdmask=tmask(:,:,jk)-ztmask0(:,:,jk); 327 DO jj = 2,jpj-1 328 DO ji = fs_2,fs_jpim1 329 jip1=ji+1; jim1=ji-1; 330 jjp1=jj+1; jjm1=jj-1; 331 summsk= (ztmask0(jip1,jj ,jk)+ztmask0(jim1,jj ,jk)+ztmask0(ji ,jjp1,jk)+ztmask0(ji ,jjm1,jk)) 332 IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 333 !! horizontal basic extrapolation 334 tsn(ji,jj,jk,1)=( zts0(jip1,jj ,jk,1)*ztmask0(jip1,jj ,jk) & 335 & +zts0(jim1,jj ,jk,1)*ztmask0(jim1,jj ,jk) & 336 & +zts0(ji ,jjp1,jk,1)*ztmask0(ji ,jjp1,jk) & 337 & +zts0(ji ,jjm1,jk,1)*ztmask0(ji ,jjm1,jk) ) / summsk 338 tsn(ji,jj,jk,2)=( zts0(jip1,jj ,jk,2)*ztmask0(jip1,jj ,jk) & 339 & +zts0(jim1,jj ,jk,2)*ztmask0(jim1,jj ,jk) & 340 & +zts0(ji ,jjp1,jk,2)*ztmask0(ji ,jjp1,jk) & 341 & +zts0(ji ,jjm1,jk,2)*ztmask0(ji ,jjm1,jk) ) / summsk 342 ztmask1(ji,jj,jk)=1 343 ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN 344 !! vertical extrapolation if horizontal extrapolation failed 345 jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1) 346 summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1)) 347 IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN 348 tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1) & 349 & +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / summsk 350 tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1) & 351 & +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / summsk 352 ztmask1(ji,jj,jk)=1._wp 353 ENDIF 354 ENDIF 355 END DO 356 END DO 357 END DO 358 ! 359 CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.) 360 ! 361 ! update 362 zts0(:,:,:,:) = tsn(:,:,:,:) 363 ztmask0 = ztmask1 364 ! 352 365 END DO 353 366 … … 395 408 ! set mbkt and mikt to 1 in thiese location 396 409 WHERE (SUM(tmask,dim=3) == 0) 397 mbkt(:,:) =1 ; mbku(:,:)=1 ; mbkv(:,:)=1398 mikt(:,:) =1 ; miku(:,:)=1 ; mikv(:,:)=1410 mbkt(:,:) = 1 ; mbku(:,:)=1 ; mbkv(:,:) = 1 411 mikt(:,:) = 1 ; miku(:,:)=1 ; mikv(:,:) = 1 399 412 END WHERE 400 413 ! -------------------------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
