Changeset 6060 for branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM
- Timestamp:
- 2015-12-16T10:25:22+01:00 (8 years ago)
- Location:
- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 1 deleted
- 11 edited
-
daymod.F90 (modified) (2 diffs)
-
dom_oce.F90 (modified) (12 diffs)
-
domain.F90 (modified) (8 diffs)
-
domcfg.F90 (modified) (1 diff)
-
domhgr.F90 (modified) (2 diffs)
-
dommsk.F90 (modified) (7 diffs)
-
domstp.F90 (modified) (2 diffs)
-
domvvl.F90 (modified) (40 diffs)
-
domzgr.F90 (modified) (49 diffs)
-
domzgr_substitute.h90 (deleted)
-
dtatsd.F90 (modified) (2 diffs)
-
istate.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
r5563 r6060 26 26 !! 27 27 !!---------------------------------------------------------------------- 28 USE dom_oce 29 USE phycst 30 USE in_out_manager 31 USE iom 32 USE ioipsl , ONLY : ymds2ju ! for calendar33 USE prtctl 34 USE trc_oce , ONLY : lk_offline ! offline flag35 USE timing 36 USE restart 28 USE dom_oce ! ocean space and time domain 29 USE phycst ! physical constants 30 USE in_out_manager ! I/O manager 31 USE iom ! 32 USE ioipsl , ONLY : ymds2ju ! for calendar 33 USE prtctl ! Print control 34 USE trc_oce , ONLY : lk_offline ! offline flag 35 USE timing ! Timing 36 USE restart ! restart 37 37 38 38 IMPLICIT NONE … … 43 43 PUBLIC day_mth ! Needed by TAM 44 44 45 INTEGER, PUBLIC :: nsecd, nsecd05, ndt, ndt05 !(PUBLIC for TAM)45 INTEGER, PUBLIC :: nsecd, nsecd05, ndt, ndt05 !: (PUBLIC for TAM) 46 46 47 47 !!---------------------------------------------------------------------- -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5930 r6060 11 11 !! to the optimization of BDY communications 12 12 !! 3.7 ! 2015-11 (G. Madec) introduce surface and scale factor ratio 13 !! - ! 2015-11 (G. Madec, A. Coward) time varying zgr by default 13 14 !!---------------------------------------------------------------------- 14 15 … … 47 48 !! Free surface parameters 48 49 !! ======================= 49 LOGICAL , PUBLIC :: ln_dynspg_exp!: Explicit free surface flag50 LOGICAL , PUBLIC :: ln_dynspg_ts!: Split-Explicit free surface flag50 LOGICAL , PUBLIC :: ln_dynspg_exp !: Explicit free surface flag 51 LOGICAL , PUBLIC :: ln_dynspg_ts !: Split-Explicit free surface flag 51 52 52 53 !! Time splitting parameters … … 61 62 !! Horizontal grid parameters for domhgr 62 63 !! ===================================== 63 INTEGER :: jphgr_msh !: type of horizontal mesh64 INTEGER :: jphgr_msh !: type of horizontal mesh 64 65 ! ! = 0 curvilinear coordinate on the sphere read in coordinate.nc 65 66 ! ! = 1 geographical mesh on the sphere with regular grid-spacing … … 68 69 ! ! = 4 Mercator grid with T/U point at the equator 69 70 70 REAL(wp) :: ppglam0 71 REAL(wp) :: ppgphi0 71 REAL(wp) :: ppglam0 !: longitude of first raw and column T-point (jphgr_msh = 1) 72 REAL(wp) :: ppgphi0 !: latitude of first raw and column T-point (jphgr_msh = 1) 72 73 ! ! used for Coriolis & Beta parameters (jphgr_msh = 2 or 3) 73 REAL(wp) :: ppe1_deg 74 REAL(wp) :: ppe2_deg 75 REAL(wp) :: ppe1_m 76 REAL(wp) :: ppe2_m 74 REAL(wp) :: ppe1_deg !: zonal grid-spacing (degrees) 75 REAL(wp) :: ppe2_deg !: meridional grid-spacing (degrees) 76 REAL(wp) :: ppe1_m !: zonal grid-spacing (degrees) 77 REAL(wp) :: ppe2_m !: meridional grid-spacing (degrees) 77 78 78 79 !! Vertical grid parameter for domzgr 79 80 !! ================================== 80 REAL(wp) :: ppsur 81 REAL(wp) :: ppa0 82 REAL(wp) :: ppa1 83 REAL(wp) :: ppkth 84 REAL(wp) :: ppacr 81 REAL(wp) :: ppsur !: ORCA r4, r2 and r05 coefficients 82 REAL(wp) :: ppa0 !: (default coefficients) 83 REAL(wp) :: ppa1 !: 84 REAL(wp) :: ppkth !: 85 REAL(wp) :: ppacr !: 85 86 ! 86 87 ! If both ppa0 ppa1 and ppsur are specified to 0, then 87 88 ! they are computed from ppdzmin, pphmax , ppkth, ppacr in dom_zgr 88 REAL(wp) :: ppdzmin 89 REAL(wp) :: pphmax 89 REAL(wp) :: ppdzmin !: Minimum vertical spacing 90 REAL(wp) :: pphmax !: Maximum depth 90 91 ! 91 LOGICAL :: ldbletanh 92 REAL(wp) :: ppa2 93 REAL(wp) :: ppkth2 94 REAL(wp) :: ppacr2 92 LOGICAL :: ldbletanh !: Use/do not use double tanf function for vertical coordinates 93 REAL(wp) :: ppa2 !: Double tanh function parameters 94 REAL(wp) :: ppkth2 !: 95 REAL(wp) :: ppacr2 !: 95 96 96 97 ! !! old non-DOCTOR names still used in the model … … 106 107 REAL(wp), PUBLIC :: rdth !: depth variation of tracer step 107 108 108 ! !!! associated variables 109 INTEGER , PUBLIC :: neuler !: restart euler forward option (0=Euler) 110 REAL(wp), PUBLIC :: atfp1 !: asselin time filter coeff. (atfp1= 1-2*atfp) 109 ! !!! associated variables 110 INTEGER , PUBLIC :: neuler !: restart euler forward option (0=Euler) 111 REAL(wp), PUBLIC :: atfp1 !: asselin time filter coeff. (atfp1= 1-2*atfp) 112 111 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: rdttra !: vertical profile of tracer time step 112 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: r2dtra !: = 2*rdttra except at nit000 (=rdttra) if neuler=0 … … 177 179 !! vertical coordinate and scale factors 178 180 !! --------------------------------------------------------------------- 179 ! !!* Namelist namzgr : vertical coordinate * 180 LOGICAL, PUBLIC :: ln_zco !: z-coordinate - full step 181 LOGICAL, PUBLIC :: ln_zps !: z-coordinate - partial step 182 LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate 183 LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF 184 185 !! All coordinates 186 !! --------------- 187 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w_0 !: depth of t-points (sum of e3w) (m) 188 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0, gdepw_0 !: analytical (time invariant) depth at t-w points (m) 189 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 , e3f_0 !: analytical (time invariant) vertical scale factors at v-f 190 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 , e3u_0 !: t-u points (m) 191 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 !: analytical (time invariant) vertical scale factors at vw 192 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 , e3uw_0 !: w-uw points (m) 193 #if defined key_vvl 194 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .TRUE. !: variable grid flag 195 196 !! All coordinates 197 !! --------------- 198 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w_n !: now depth of T-points (sum of e3w) (m) 199 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_n, gdepw_n !: now depth at T-W points (m) 200 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_b, gdepw_b !: before depth at T-W points (m) 201 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_n !: now vertical scale factors at t point (m) 202 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_n , e3v_n !: - - - - u --v points (m) 203 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_n , e3f_n !: - - - - w --f points (m) 204 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_n , e3vw_n !: - - - - uw--vw points (m) 205 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_b !: before - - - - t points (m) 206 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_b !: before - - - - t points (m) 207 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_b , e3v_b !: - - - - - u --v points (m) 208 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_b , e3vw_b !: - - - - - uw--vw points (m) 209 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_a !: after - - - - t point (m) 210 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_a , e3v_a !: - - - - - u --v points (m) 211 #else 212 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .FALSE. !: fixed grid flag 213 #endif 214 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hur , hvr !: Now inverse of u and v-points ocean depth (1/m) 215 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu , hv !: depth at u- and v-points (meters) 216 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht !: depth at t-points (meters) 217 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ehur_a, ehvr_a !: After inverse of u and v-points ocean depth (1/m) 218 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ehu_a , ehv_a !: depth at u- and v-points (meters) 219 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ehur_b, ehvr_b !: Before inverse of u and v-points ocean depth (1/m) 220 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ehu_b , ehv_b !: depth at u- and v-points (meters) 221 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 !: reference depth at t- points (meters) 222 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hv_0 !: reference depth at u- and v-points (meters) 181 ! !!* Namelist namzgr : vertical coordinate * 182 LOGICAL, PUBLIC :: ln_zco !: z-coordinate - full step 183 LOGICAL, PUBLIC :: ln_zps !: z-coordinate - partial step 184 LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate 185 LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF 186 LOGICAL, PUBLIC :: ln_linssh !: variable grid flag 187 188 ! ! ref. ! before ! now ! after ! 189 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 , e3t_b , e3t_n , e3t_a !: t- vert. scale factor [m] 190 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_0 , e3u_b , e3u_n , e3u_a !: u- vert. scale factor [m] 191 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 , e3v_b , e3v_n , e3v_a !: v- vert. scale factor [m] 192 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f_0 , e3f_n !: f- vert. scale factor [m] 193 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 , e3w_b , e3w_n !: w- vert. scale factor [m] 194 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_0 , e3uw_b , e3uw_n !: uw-vert. scale factor [m] 195 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 , e3vw_b , e3vw_n !: vw-vert. scale factor [m] 196 197 ! ! ref. ! before ! now ! 198 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 , gdept_b , gdept_n !: t- depth [m] 199 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdepw_0 , gdepw_b , gdepw_n !: w- depth [m] 200 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0 , gde3w_n !: w- depth (sum of e3w) [m] 201 202 ! ! ref. ! before ! now ! after ! 203 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 , ht_n !: t-depth [m] 204 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hu_b , hu_n , hu_a !: u-depth [m] 205 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0 , hv_b , hv_n , hv_a !: u-depth [m] 206 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_b , r1_hu_n , r1_hu_a !: inverse of u-depth [1/m] 207 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_b , r1_hv_n , r1_hv_a !: inverse of v-depth [1/m] 208 223 209 224 210 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) 225 211 INTEGER, PUBLIC :: nlb10 !: shallowest W level Bellow ~10m (nla10 + 1) 226 212 227 !! z-coordinate with full steps (also used in the other cases as reference z-coordinate)213 !! 1D reference vertical coordinate 228 214 !! =-----------------====------ 229 215 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m) … … 231 217 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3tp , e3wp !: ocean bottom level thickness at T and W points 232 218 219 !!gm This should be removed from here.... ==>>> only used in domzgr at initialization phase 233 220 !! s-coordinate and hybrid z-s-coordinate 234 221 !! =----------------======--------------- … … 244 231 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing t--u points (m) 245 232 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio 233 !!gm end 246 234 247 235 !!---------------------------------------------------------------------- … … 249 237 !! --------------------------------------------------------------------- 250 238 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbathy !: number of ocean level (=0, 1, ... , jpk-1) 251 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level 252 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level 253 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 239 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt, mbku, mbkv !: vertical index of the bottom last T-, U- & V ocean level 240 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 254 241 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask 255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function256 242 257 243 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfdep !: top first ocean level (ISF) … … 352 338 & ff (jpi,jpj) , STAT=ierr(3) ) 353 339 ! 354 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) , & 355 & gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , & 356 & gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 357 ! 358 #if defined key_vvl 359 ALLOCATE( gdep3w_n(jpi,jpj,jpk) , e3t_n (jpi,jpj,jpk) , e3u_n (jpi,jpj,jpk) , & 360 & gdept_n (jpi,jpj,jpk) , e3v_n (jpi,jpj,jpk) , e3w_n (jpi,jpj,jpk) , & 361 & gdepw_n (jpi,jpj,jpk) , e3f_n (jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , e3uw_n(jpi,jpj,jpk) , & 362 & e3t_b (jpi,jpj,jpk) , e3u_b (jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) , & 363 & e3uw_b (jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) , & 364 & gdept_b (jpi,jpj,jpk) ,gdepw_b(jpi,jpj,jpk) , e3w_b (jpi,jpj,jpk) , & 365 & e3t_a (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk) , & 366 & ehu_a (jpi,jpj) , ehv_a (jpi,jpj), & 367 & ehur_a (jpi,jpj) , ehvr_a(jpi,jpj), & 368 & ehu_b (jpi,jpj) , ehv_b (jpi,jpj), & 369 & ehur_b (jpi,jpj) , ehvr_b(jpi,jpj), STAT=ierr(5) ) 370 #endif 371 ! 372 ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) , ht_0(jpi,jpj) , & 373 & hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , ht (jpi,jpj) , STAT=ierr(6) ) 340 ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) , & 341 & gdept_b(jpi,jpj,jpk) , gdepw_b(jpi,jpj,jpk) , & 342 & gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) ) 343 ! 344 ALLOCATE( e3t_0(jpi,jpj,jpk) , e3u_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , & 345 & e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) , e3w_b(jpi,jpj,jpk) , & & 346 & e3t_n(jpi,jpj,jpk) , e3u_n(jpi,jpj,jpk) , e3v_n(jpi,jpj,jpk) , e3f_n(jpi,jpj,jpk) , e3w_n(jpi,jpj,jpk) , & & 347 & e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) , & 348 ! ! 349 & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , & 350 & e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) , & 351 & e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , STAT=ierr(5) ) 352 ! 353 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , & 354 & hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) , & 355 & ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) , & 356 & hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6) ) 357 ! 374 358 ! 375 359 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , & … … 386 370 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 387 371 & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 388 & bmask (jpi,jpj) , &389 372 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 390 373 … … 405 388 !!====================================================================== 406 389 END MODULE dom_oce 407 -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5836 r6060 13 13 !! 3.3 ! 2010-11 (G. Madec) initialisation in C1D configuration 14 14 !! 3.6 ! 2013 ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs 15 !! 3.7 ! 2015-11 (G. Madec, A. Coward) time varying zgr by default 15 16 !!---------------------------------------------------------------------- 16 17 … … 36 37 ! 37 38 USE in_out_manager ! I/O manager 39 USE wrk_nemo ! Memory Allocation 38 40 USE lib_mpp ! distributed memory computing library 39 41 USE lbclnk ! ocean lateral boundary condition (or mpp link) … … 45 47 PUBLIC dom_init ! called by opa.F90 46 48 47 !! * Substitutions48 # include "domzgr_substitute.h90"49 49 !!------------------------------------------------------------------------- 50 50 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 70 70 !! - 1D configuration, move Coriolis, u and v at T-point 71 71 !!---------------------------------------------------------------------- 72 INTEGER :: jk ! dummy loop argument72 INTEGER :: jk ! dummy loop indices 73 73 INTEGER :: iconf = 0 ! local integers 74 REAL(wp), POINTER, DIMENSION(:,:) :: z1_hu_0, z1_hv_0 74 75 !!---------------------------------------------------------------------- 75 76 ! … … 82 83 ENDIF 83 84 ! 84 CALL dom_nam ! read namelist ( namrun, namdom ) 85 CALL dom_clo ! Closed seas and lake 86 CALL dom_hgr ! Horizontal mesh 87 CALL dom_zgr ! Vertical mesh and bathymetry 88 CALL dom_msk ! Masks 89 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 90 ! 91 ht_0(:,:) = 0._wp ! Reference ocean depth at T-points 92 hu_0(:,:) = 0._wp ! Reference ocean depth at U-points 93 hv_0(:,:) = 0._wp ! Reference ocean depth at V-points 94 DO jk = 1, jpk 85 ! !== Reference coordinate system ==! 86 ! 87 CALL dom_nam ! read namelist ( namrun, namdom ) 88 CALL dom_clo ! Closed seas and lake 89 CALL dom_hgr ! Horizontal mesh 90 CALL dom_zgr ! Vertical mesh and bathymetry 91 CALL dom_msk ! Masks 92 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 93 ! 94 ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1) ! Reference ocean thickness 95 hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1) 96 hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1) 97 DO jk = 2, jpk 95 98 ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 96 99 hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) … … 98 101 END DO 99 102 ! 100 IF( lk_vvl ) CALL dom_vvl_init ! Vertical variable mesh 101 ! 102 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 103 ! 104 ! 105 hu(:,:) = 0._wp ! Ocean depth at U-points 106 hv(:,:) = 0._wp ! Ocean depth at V-points 107 ht(:,:) = 0._wp ! Ocean depth at T-points 108 DO jk = 1, jpkm1 109 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 110 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 111 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 112 END DO 113 ! ! Inverse of the local depth 114 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 115 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 116 117 CALL dom_stp ! time step 118 IF( nmsh /= 0 ) CALL dom_wri ! Create a domain file 119 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 103 ! !== time varying part of coordinate system ==! 104 ! 105 IF( ln_linssh ) THEN ! Fix in time : set to the reference one for all 106 ! before ! now ! after ! 107 ; gdept_b = gdept_0 ; gdept_n = gdept_0 ! --- ! depth of grid-points 108 ; gdepw_b = gdepw_0 ; gdepw_n = gdepw_0 ! --- ! 109 ; ; gde3w_n = gde3w_0 ! --- ! 110 ! 111 ; e3t_b = e3t_0 ; e3t_n = e3t_0 ; e3t_a = e3t_0 ! scale factors 112 ; e3u_b = e3u_0 ; e3u_n = e3u_0 ; e3u_a = e3u_0 ! 113 ; e3v_b = e3v_0 ; e3v_n = e3v_0 ; e3v_a = e3v_0 ! 114 ; ; e3f_n = e3f_0 ! --- ! 115 ; e3w_b = e3w_0 ; e3w_n = e3w_0 ! --- ! 116 ; e3uw_b = e3uw_0 ; e3uw_n = e3uw_0 ! --- ! 117 ; e3vw_b = e3vw_0 ; e3vw_n = e3vw_0 ! --- ! 118 ! 119 CALL wrk_alloc( jpi,jpj, z1_hu_0, z1_hv_0 ) 120 ! 121 z1_hu_0(:,:) = umask_i(:,:) / ( hu_0(:,:) + 1._wp - umask_i(:,:) ) ! _i mask due to ISF 122 z1_hv_0(:,:) = vmask_i(:,:) / ( hv_0(:,:) + 1._wp - vmask_i(:,:) ) 123 ! 124 ! before ! now ! after ! 125 ; ; ht_n = ht_0 ! ! water column thickness 126 ; hu_b = hu_0 ; hu_n = hu_0 ; hu_a = hu_0 ! 127 ; hv_b = hv_0 ; hv_n = hv_0 ; hv_a = hv_0 ! 128 ; r1_hu_b = z1_hu_0 ; r1_hu_n = z1_hu_0 ; r1_hu_a = z1_hu_0 ! inverse of water column thickness 129 ; r1_hv_b = z1_hv_0 ; r1_hv_n = z1_hv_0 ; r1_hv_a = z1_hv_0 ! 130 ! 131 CALL wrk_dealloc( jpi,jpj, z1_hu_0, z1_hv_0 ) 132 ! 133 ELSE ! time varying : initialize before/now/after variables 134 ! 135 CALL dom_vvl_init 136 ! 137 ENDIF 138 ! 139 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 140 ! 141 CALL dom_stp ! time step 142 IF( nmsh /= 0 ) CALL dom_wri ! Create a domain file 143 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 120 144 ! 121 145 IF( nn_timing == 1 ) CALL timing_stop('dom_init') … … 403 427 INTEGER :: ji, jj, jk 404 428 REAL(wp) :: zrxmax 405 REAL(wp), DIMENSION(4) :: zr1429 REAL(wp), DIMENSION(4) :: zr1 406 430 !!---------------------------------------------------------------------- 407 431 rx1(:,:) = 0._wp … … 412 436 DO jj = 2, jpjm1 413 437 DO jk = 1, jpkm1 414 zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw_0(ji ,jj ,jk )-gdepw_0(ji-1,jj ,jk )&415 & +gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji-1,jj ,jk+1))&416 & /(gdepw_0(ji ,jj ,jk )+gdepw_0(ji-1,jj ,jk )&417 & -gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji-1,jj ,jk+1) + rsmall))418 zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw_0(ji+1,jj ,jk )-gdepw_0(ji ,jj ,jk )&419 & +gdepw_0(ji+1,jj ,jk+1)-gdepw_0(ji ,jj ,jk+1))&420 & /(gdepw_0(ji+1,jj ,jk )+gdepw_0(ji ,jj ,jk )&421 & -gdepw_0(ji+1,jj ,jk+1)-gdepw_0(ji ,jj ,jk+1) + rsmall))422 zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw_0(ji ,jj+1,jk )-gdepw_0(ji ,jj ,jk )&423 & +gdepw_0(ji ,jj+1,jk+1)-gdepw_0(ji ,jj ,jk+1))&424 & /(gdepw_0(ji ,jj+1,jk )+gdepw_0(ji ,jj ,jk )&425 & -gdepw_0(ji ,jj+1,jk+1)-gdepw_0(ji ,jj ,jk+1) + rsmall))426 zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw_0(ji ,jj ,jk )-gdepw_0(ji ,jj-1,jk )&427 & +gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji ,jj-1,jk+1))&428 & /(gdepw_0(ji ,jj ,jk )+gdepw_0(ji ,jj-1,jk )&429 & -gdepw_0(ji, jj ,jk+1)-gdepw_0(ji ,jj-1,jk+1) + rsmall))430 zrxmax = MAXVAL( zr1(1:4))431 rx1(ji,jj) = MAX( rx1(ji,jj), zrxmax)438 zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) & 439 & +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) & 440 & / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) & 441 & -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk) 442 zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) & 443 & +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) & 444 & / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) & 445 & -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk) 446 zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) & 447 & +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) & 448 & / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) & 449 & -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk) 450 zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) & 451 & +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) & 452 & / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) & 453 & -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk) 454 zrxmax = MAXVAL( zr1(1:4) ) 455 rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax ) 432 456 END DO 433 457 END DO -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90
r4667 r6060 118 118 WRITE(numout,*) ' south-west indices jpizoom = ', jpizoom, & 119 119 & ' jpjzoom = ', jpjzoom 120 WRITE(numout,*) 121 WRITE(numout,*) ' conversion local ==> data i-index domain' 122 WRITE(numout,25) (mig(ji),ji = 1,jpi) 123 WRITE(numout,*) 124 WRITE(numout,*) ' conversion data ==> local i-index domain' 125 WRITE(numout,*) ' starting index' 126 WRITE(numout,25) (mi0(ji),ji = 1,jpidta) 127 WRITE(numout,*) ' ending index' 128 WRITE(numout,25) (mi1(ji),ji = 1,jpidta) 129 WRITE(numout,*) 130 WRITE(numout,*) ' conversion local ==> data j-index domain' 131 WRITE(numout,25) (mjg(jj),jj = 1,jpj) 132 WRITE(numout,*) 133 WRITE(numout,*) ' conversion data ==> local j-index domain' 134 WRITE(numout,*) ' starting index' 135 WRITE(numout,25) (mj0(jj),jj = 1,jpjdta) 136 WRITE(numout,*) ' ending index' 137 WRITE(numout,25) (mj1(jj),jj = 1,jpjdta) 120 IF( nn_print >= 1 ) THEN 121 WRITE(numout,*) 122 WRITE(numout,*) ' conversion local ==> data i-index domain' 123 WRITE(numout,25) (mig(ji),ji = 1,jpi) 124 WRITE(numout,*) 125 WRITE(numout,*) ' conversion data ==> local i-index domain' 126 WRITE(numout,*) ' starting index' 127 WRITE(numout,25) (mi0(ji),ji = 1,jpidta) 128 WRITE(numout,*) ' ending index' 129 WRITE(numout,25) (mi1(ji),ji = 1,jpidta) 130 WRITE(numout,*) 131 WRITE(numout,*) ' conversion local ==> data j-index domain' 132 WRITE(numout,25) (mjg(jj),jj = 1,jpj) 133 WRITE(numout,*) 134 WRITE(numout,*) ' conversion data ==> local j-index domain' 135 WRITE(numout,*) ' starting index' 136 WRITE(numout,25) (mj0(jj),jj = 1,jpjdta) 137 WRITE(numout,*) ' ending index' 138 WRITE(numout,25) (mj1(jj),jj = 1,jpjdta) 139 ENDIF 138 140 ENDIF 139 141 25 FORMAT( 100(10x,19i4,/) ) -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r5836 r6060 348 348 e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 349 349 350 IF( lwp .AND. .NOT.ln_rstart ) THEN ! Control print : Grid informations (if not restart)350 IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN ! Control print : Grid informations (if not restart) 351 351 WRITE(numout,*) 352 352 WRITE(numout,*) ' longitude and e1 scale factors' … … 391 391 ! 392 392 #if defined key_agrif 393 IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN! for EEL6 configuration only394 IF( .NOT. 393 IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only 394 IF( .NOT.Agrif_Root() ) THEN 395 395 zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 396 396 ENDIF -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r5930 r6060 7 7 !! 6.0 ! 1993-03 (M. Guyon) symetrical conditions (M. Guyon) 8 8 !! 7.0 ! 1996-01 (G. Madec) suppression of common work arrays 9 !! - ! 1996-05 (G. Madec) mask computed from tmask and sup- 10 !! ! pression of the double computation of bmask 9 !! - ! 1996-05 (G. Madec) mask computed from tmask 11 10 !! 8.0 ! 1997-02 (G. Madec) mesh information put in domhgr.F 12 11 !! 8.1 ! 1997-07 (G. Madec) modification of mbathy and fmask … … 25 24 USE oce ! ocean dynamics and tracers 26 25 USE dom_oce ! ocean space and time domain 26 ! 27 27 USE in_out_manager ! I/O manager 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 USE lib_mpp 29 USE lib_mpp ! 30 30 USE wrk_nemo ! Memory allocation 31 31 USE timing ! Timing … … 34 34 PRIVATE 35 35 36 PUBLIC dom_msk 36 PUBLIC dom_msk ! routine called by inidom.F90 37 37 38 38 ! !!* Namelist namlbc : lateral boundary condition * … … 89 89 !! 90 90 !! N.B. If nperio not equal to 0, the land/ocean mask arrays 91 !! are defined with the proper value at lateral domain boundaries, 92 !! but bmask. indeed, bmask defined the domain over which the 93 !! barotropic stream function is computed. this domain cannot 94 !! contain identical columns because the matrix associated with 95 !! the barotropic stream function equation is then no more inverti- 96 !! ble. therefore bmask is set to 0 along lateral domain boundaries 97 !! even IF nperio is not zero. 91 !! are defined with the proper value at lateral domain boundaries. 98 92 !! 99 93 !! In case of open boundaries (lk_bdy=T): 100 94 !! - tmask is set to 1 on the points to be computed bay the open 101 95 !! boundaries routines. 102 !! - bmask is set to 0 on the open boundaries.103 96 !! 104 97 !! ** Action : tmask : land/ocean mask at t-point (=0. or 1.) … … 107 100 !! fmask : land/ocean mask at f-point (=0. or 1.) 108 101 !! =rn_shlat along lateral boundaries 109 !! bmask : land/ocean mask at barotropic stream110 !! function point (=0. or 1.) and set to 0 along lateral boundaries111 102 !! tmask_i : interior ocean mask 112 103 !!---------------------------------------------------------------------- … … 254 245 END DO 255 246 256 ! 4. ocean/land mask for the elliptic equation257 ! --------------------------------------------258 bmask(:,:) = ssmask(:,:) ! elliptic equation is written at t-point259 !260 ! ! Boundary conditions261 ! ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi262 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN263 bmask( 1 ,:) = 0._wp264 bmask(jpi,:) = 0._wp265 ENDIF266 IF( nperio == 2 ) THEN ! south symmetric : bmask must be set to 0. on row 1267 bmask(:, 1 ) = 0._wp268 ENDIF269 ! ! north fold :270 IF( nperio == 3 .OR. nperio == 4 ) THEN ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row271 DO ji = 1, jpi272 ii = ji + nimpp - 1273 bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)274 bmask(ji,jpj ) = 0._wp275 END DO276 ENDIF277 IF( nperio == 5 .OR. nperio == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj278 bmask(:,jpj) = 0._wp279 ENDIF280 !281 IF( lk_mpp ) THEN ! mpp specificities282 ! ! bmask is set to zero on the overlap region283 IF( nbondi /= -1 .AND. nbondi /= 2 ) bmask( 1 :jpreci,:) = 0._wp284 IF( nbondi /= 1 .AND. nbondi /= 2 ) bmask(nlci:jpi ,:) = 0._wp285 IF( nbondj /= -1 .AND. nbondj /= 2 ) bmask(:, 1 :jprecj) = 0._wp286 IF( nbondj /= 1 .AND. nbondj /= 2 ) bmask(:,nlcj:jpj ) = 0._wp287 !288 IF( npolj == 3 .OR. npolj == 4 ) THEN ! north fold : bmask must be set to 0. on rows jpj-1 and jpj289 DO ji = 1, nlci290 ii = ji + nimpp - 1291 bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)292 bmask(ji,nlcj ) = 0._wp293 END DO294 ENDIF295 IF( npolj == 5 .OR. npolj == 6 ) THEN ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj296 DO ji = 1, nlci297 bmask(ji,nlcj ) = 0._wp298 END DO299 ENDIF300 ENDIF301 302 247 ! Lateral boundary conditions on velocity (modify fmask) 303 248 ! --------------------------------------- … … 399 344 ! 400 345 CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask 401 346 ! 402 347 ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) 403 404 IF( nprint == 1 .AND. lwp ) THEN ! Control print405 imsk(:,:) = INT( tmask_i(:,:) )406 WRITE(numout,*) ' tmask_i : '407 CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &408 & 1, jpj, 1, 1, numout)409 WRITE (numout,*)410 WRITE (numout,*) ' dommsk: tmask for each level'411 WRITE (numout,*) ' ----------------------------'412 DO jk = 1, jpk413 imsk(:,:) = INT( tmask(:,:,jk) )414 415 WRITE(numout,*)416 WRITE(numout,*) ' level = ',jk417 CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &418 & 1, jpj, 1, 1, numout)419 END DO420 WRITE(numout,*)421 WRITE(numout,*) ' dom_msk: vmask for each level'422 WRITE(numout,*) ' -----------------------------'423 DO jk = 1, jpk424 imsk(:,:) = INT( vmask(:,:,jk) )425 WRITE(numout,*)426 WRITE(numout,*) ' level = ',jk427 CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &428 & 1, jpj, 1, 1, numout)429 END DO430 WRITE(numout,*)431 WRITE(numout,*) ' dom_msk: fmask for each level'432 WRITE(numout,*) ' -----------------------------'433 DO jk = 1, jpk434 imsk(:,:) = INT( fmask(:,:,jk) )435 WRITE(numout,*)436 WRITE(numout,*) ' level = ',jk437 CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &438 & 1, jpj, 1, 1, numout )439 END DO440 WRITE(numout,*)441 WRITE(numout,*) ' dom_msk: bmask '442 WRITE(numout,*) ' ---------------'443 WRITE(numout,*)444 imsk(:,:) = INT( bmask(:,:) )445 CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, &446 & 1, jpj, 1, 1, numout )447 ENDIF448 348 ! 449 349 CALL wrk_dealloc( jpi, jpj, imsk ) -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domstp.F90
r4292 r6060 22 22 PUBLIC dom_stp ! routine called by inidom.F90 23 23 24 !! * Substitutions25 # include "domzgr_substitute.h90"26 24 !!---------------------------------------------------------------------- 27 25 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 84 82 IF(lwp) WRITE(numout,*)' accelerating the convergence' 85 83 IF(lwp) WRITE(numout,*)' dynamics time step = ', rdt/3600., ' hours' 86 IF( ln_sco .AND. rdtmin /= rdtmax .AND. lk_vvl) &84 IF( ln_sco .AND. rdtmin /= rdtmax .AND. .NOT.ln_linssh ) & 87 85 & CALL ctl_stop ( ' depth dependent acceleration of convergence not implemented in s-coordinates & 88 86 & nor in variable volume' ) -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5836 r6060 20 20 !!---------------------------------------------------------------------- 21 21 USE oce ! ocean dynamics and tracers 22 USE phycst ! physical constant 22 23 USE dom_oce ! ocean space and time domain 23 24 USE sbc_oce ! ocean surface boundary condition 25 USE restart ! ocean restart 26 ! 24 27 USE in_out_manager ! I/O manager 25 28 USE iom ! I/O manager library 26 USE restart ! ocean restart27 29 USE lib_mpp ! distributed memory computing library 28 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 60 62 61 63 !! * Substitutions 62 # include "domzgr_substitute.h90"63 64 # include "vectopt_loop_substitute.h90" 64 65 !!---------------------------------------------------------------------- 65 !! NEMO/OPA 3. 3 , NEMO-Consortium (2010)66 !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 66 67 !! $Id$ 67 68 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 68 69 !!---------------------------------------------------------------------- 69 70 70 CONTAINS 71 71 … … 74 74 !! *** FUNCTION dom_vvl_alloc *** 75 75 !!---------------------------------------------------------------------- 76 IF( ln_vvl_zstar ) dom_vvl_alloc = 076 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 77 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 78 78 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & … … 81 81 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 82 82 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 83 un_td = 0. 0_wp84 vn_td = 0. 0_wp83 un_td = 0._wp 84 vn_td = 0._wp 85 85 ENDIF 86 86 IF( ln_vvl_ztilde ) THEN … … 89 89 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 90 90 ENDIF 91 91 ! 92 92 END FUNCTION dom_vvl_alloc 93 93 … … 103 103 !! - interpolate scale factors 104 104 !! 105 !! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b)106 !! - Regrid: fse3(u/v)_n107 !! fse3(u/v)_b108 !! fse3w_n109 !! fse3(u/v)w_b110 !! fse3(u/v)w_n111 !! fsdept_n, fsdepw_n and fsde3w_n105 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 106 !! - Regrid: e3(u/v)_n 107 !! e3(u/v)_b 108 !! e3w_n 109 !! e3(u/v)w_b 110 !! e3(u/v)w_n 111 !! gdept_n, gdepw_n and gde3w_n 112 112 !! - h(t/u/v)_0 113 113 !! - frq_rst_e3t and frq_rst_hdv … … 115 115 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 116 !!---------------------------------------------------------------------- 117 USE phycst, ONLY : rpi, rsmall, rad 118 !! * Local declarations 119 INTEGER :: ji,jj,jk 117 INTEGER :: ji, jj, jk 120 118 INTEGER :: ii0, ii1, ij0, ij1 121 119 REAL(wp):: zcoef 122 120 !!---------------------------------------------------------------------- 123 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 124 121 ! 122 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 123 ! 125 124 IF(lwp) WRITE(numout,*) 126 125 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 127 126 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 128 129 ! choose vertical coordinate (z_star, z_tilde or layer) 130 ! ========================== 131 CALL dom_vvl_ctl 132 133 ! Allocate module arrays 134 ! ====================== 127 ! 128 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 129 ! 130 ! ! Allocate module arrays 135 131 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 136 137 ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk)) 138 ! ============================================================================ 132 ! 133 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 139 134 CALL dom_vvl_rst( nit000, 'READ' ) 140 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 141 142 ! Reconstruction of all vertical scale factors at now and before time steps 143 ! ============================================================================= 144 ! Horizontal scale factor interpolations 145 ! -------------------------------------- 146 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 147 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 148 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 149 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 150 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 151 ! Vertical scale factor interpolations 152 ! ------------------------------------ 153 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 154 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 155 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 156 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 157 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 158 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 159 ! t- and w- points depth 160 ! ---------------------- 161 ! set the isf depth as it is in the initial step 162 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 163 fsdepw_n(:,:,1) = 0.0_wp 164 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 165 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 166 fsdepw_b(:,:,1) = 0.0_wp 167 168 DO jk = 2, jpk 135 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 136 ! 137 ! !== Set of all other vertical scale factors ==! (now and before) 138 ! ! Horizontal interpolation of e3t 139 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U 140 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 141 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V 142 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 143 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F 144 ! ! Vertical interpolation of e3t,u,v 145 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W 146 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' ) 147 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW 148 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 149 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW 150 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 151 ! 152 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 153 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 154 gdepw_n(:,:,1) = 0.0_wp 155 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 156 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 157 gdepw_b(:,:,1) = 0.0_wp 158 DO jk = 2, jpk ! vertical sum 169 159 DO jj = 1,jpj 170 160 DO ji = 1,jpi 171 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 172 ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 173 ! 0.5 where jk = mikt 174 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 175 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 176 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 177 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 178 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 179 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 180 fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) & 181 & + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)) 161 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 162 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 163 ! ! 0.5 where jk = mikt 164 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 165 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 166 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 167 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 168 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 169 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 170 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 171 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 172 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 182 173 END DO 183 174 END DO 184 175 END DO 185 186 ! Before depth and Inverse of the local depth of the water column at u- and v- points 187 ! ----------------------------------------------------------------------------------- 188 hu_b(:,:) = 0. 189 hv_b(:,:) = 0. 190 DO jk = 1, jpkm1 191 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 192 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 176 ! 177 ! !== thickness of the water column !! (ocean portion only) 178 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 179 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 180 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 181 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 182 hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1) 183 DO jk = 2, jpkm1 184 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 185 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 186 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 187 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 188 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 193 189 END DO 194 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 195 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 196 197 ! Restoring frequencies for z_tilde coordinate 198 ! ============================================ 190 ! 191 ! !== inverse of water column thickness ==! (u- and v- points) 192 r1_hu_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) ! _i mask due to ISF 193 r1_hu_n(:,:) = umask_i(:,:) / ( hu_n(:,:) + 1._wp - umask_i(:,:) ) 194 r1_hv_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 195 r1_hv_n(:,:) = vmask_i(:,:) / ( hv_n(:,:) + 1._wp - vmask_i(:,:) ) 196 197 ! !== z_tilde coordinate case ==! (Restoring frequencies) 199 198 IF( ln_vvl_ztilde ) THEN 200 ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 201 frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 202 frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 203 IF( ln_vvl_ztilde_as_zstar ) THEN 204 ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 205 frq_rst_e3t(:,:) = 0.0_wp 206 frq_rst_hdv(:,:) = 1.0_wp / rdt 207 ENDIF 208 IF ( ln_vvl_zstar_at_eqtor ) THEN 199 !!gm : idea: add here a READ in a file of custumized restoring frequency 200 ! ! Values in days provided via the namelist 201 ! ! use rsmall to avoid possible division by zero errors with faulty settings 202 frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 203 frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 204 ! 205 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 206 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 207 frq_rst_hdv(:,:) = 1._wp / rdt 208 ENDIF 209 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 209 210 DO jj = 1, jpj 210 211 DO ji = 1, jpi 212 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 211 213 IF( ABS(gphit(ji,jj)) >= 6.) THEN 212 214 ! values outside the equatorial band and transition zone (ztilde) 213 215 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 214 216 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 215 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 217 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 216 218 ! values inside the equatorial band (ztilde as zstar) 217 219 frq_rst_e3t(ji,jj) = 0.0_wp 218 220 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 219 ELSE 220 ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)221 ELSE ! transition band (2.5 to 6 degrees N/S) 222 ! ! (linearly transition from z-tilde to z-star) 221 223 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 222 224 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & … … 229 231 END DO 230 232 END DO 231 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN 232 ii0 = 103 ; ii1 = 111 ! Suppress ztilde in the Foxe Basin for ORCA2233 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 234 ii0 = 103 ; ii1 = 111 233 235 ij0 = 128 ; ij1 = 135 ; 234 236 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp … … 237 239 ENDIF 238 240 ENDIF 239 241 ! 240 242 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 241 243 ! 242 244 END SUBROUTINE dom_vvl_init 243 245 … … 261 263 !! - tilde_e3t_a: after increment of vertical scale factor 262 264 !! in z_tilde case 263 !! - fse3(t/u/v)_a265 !! - e3(t/u/v)_a 264 266 !! 265 267 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 266 268 !!---------------------------------------------------------------------- 267 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 268 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 269 !! * Arguments 270 INTEGER, INTENT( in ) :: kt ! time step 271 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 272 !! * Local declarations 273 INTEGER :: ji, jj, jk ! dummy loop indices 274 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 275 REAL(wp) :: z2dt ! temporary scalars 276 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 277 LOGICAL :: ll_do_bclinic ! temporary logical 278 !!---------------------------------------------------------------------- 279 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 280 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 281 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 282 283 IF(kt == nit000) THEN 269 INTEGER, INTENT( in ) :: kt ! time step 270 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 271 ! 272 INTEGER :: ji, jj, jk ! dummy loop indices 273 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 274 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars 275 LOGICAL :: ll_do_bclinic ! local logical 276 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 277 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 278 !!---------------------------------------------------------------------- 279 ! 280 IF( ln_linssh ) RETURN ! No calculation in linear free surface 281 ! 282 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 283 ! 284 CALL wrk_alloc( jpi,jpj,zht, z_scale, zwu, zwv, zhdiv ) 285 CALL wrk_alloc( jpi,jpj,jpk, ze3t ) 286 287 IF( kt == nit000 ) THEN 284 288 IF(lwp) WRITE(numout,*) 285 289 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' … … 289 293 ll_do_bclinic = .TRUE. 290 294 IF( PRESENT(kcall) ) THEN 291 IF ( kcall == 2 .AND. ln_vvl_ztilde )ll_do_bclinic = .FALSE.295 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 292 296 ENDIF 293 297 … … 295 299 ! After acale factors at t-points ! 296 300 ! ******************************* ! 297 298 301 ! ! --------------------------------------------- ! 299 302 ! ! z_star coordinate and barotropic z-tilde part ! 300 303 ! ! --------------------------------------------- ! 301 304 ! 302 305 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 303 306 DO jk = 1, jpkm1 304 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)305 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)307 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 308 e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 306 309 END DO 307 310 ! 308 311 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 309 312 ! ! ------baroclinic part------ ! 310 311 313 ! I - initialization 312 314 ! ================== … … 314 316 ! 1 - barotropic divergence 315 317 ! ------------------------- 316 zhdiv(:,:) = 0. 317 zht(:,:) = 0. 318 zhdiv(:,:) = 0._wp 319 zht(:,:) = 0._wp 318 320 DO jk = 1, jpkm1 319 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)320 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)321 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 322 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 321 323 END DO 322 324 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 325 327 ! -------------------------------------------------- 326 328 IF( ln_vvl_ztilde ) THEN 327 IF( kt .GT.nit000 ) THEN329 IF( kt > nit000 ) THEN 328 330 DO jk = 1, jpkm1 329 331 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 330 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )332 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 331 333 END DO 332 334 ENDIF 333 END 335 ENDIF 334 336 335 337 ! II - after z_tilde increments of vertical scale factors 336 338 ! ======================================================= 337 tilde_e3t_a(:,:,:) = 0. 0_wp ! tilde_e3t_a used to store tendency terms339 tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms 338 340 339 341 ! 1 - High frequency divergence term … … 341 343 IF( ln_vvl_ztilde ) THEN ! z_tilde case 342 344 DO jk = 1, jpkm1 343 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )345 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 344 346 END DO 345 347 ELSE ! layer case 346 348 DO jk = 1, jpkm1 347 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)348 END DO 349 END 349 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 350 END DO 351 ENDIF 350 352 351 353 ! 2 - Restoring term (z-tilde case only) … … 355 357 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 356 358 END DO 357 END 359 ENDIF 358 360 359 361 ! 3 - Thickness diffusion term 360 362 ! ---------------------------- 361 zwu(:,:) = 0.0_wp 362 zwv(:,:) = 0.0_wp 363 ! a - first derivative: diffusive fluxes 364 DO jk = 1, jpkm1 363 zwu(:,:) = 0._wp 364 zwv(:,:) = 0._wp 365 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 365 366 DO jj = 1, jpjm1 366 367 DO ji = 1, fs_jpim1 ! vector opt. … … 374 375 END DO 375 376 END DO 376 ! b - correction for last oceanic u-v points 377 DO jj = 1, jpj 377 DO jj = 1, jpj ! b - correction for last oceanic u-v points 378 378 DO ji = 1, jpi 379 379 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) … … 381 381 END DO 382 382 END DO 383 ! c - second derivative: divergence of diffusive fluxes 384 DO jk = 1, jpkm1 383 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 385 384 DO jj = 2, jpjm1 386 385 DO ji = fs_2, fs_jpim1 ! vector opt. … … 391 390 END DO 392 391 END DO 393 ! d - thickness diffusion transport: boundary conditions394 ! (stored for tracer advction and continuity equation)392 ! ! d - thickness diffusion transport: boundary conditions 393 ! (stored for tracer advction and continuity equation) 395 394 CALL lbc_lnk( un_td , 'U' , -1._wp) 396 395 CALL lbc_lnk( vn_td , 'V' , -1._wp) … … 410 409 ! Maximum deformation control 411 410 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 412 ze3t(:,:,jpk) = 0. 0_wp411 ze3t(:,:,jpk) = 0._wp 413 412 DO jk = 1, jpkm1 414 413 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) … … 462 461 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 463 462 DO jk = 1, jpkm1 464 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)463 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 465 464 END DO 466 465 … … 470 469 ! ! ---baroclinic part--------- ! 471 470 DO jk = 1, jpkm1 472 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)471 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 473 472 END DO 474 473 ENDIF … … 485 484 zht(:,:) = 0.0_wp 486 485 DO jk = 1, jpkm1 487 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)486 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 488 487 END DO 489 488 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 490 489 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 491 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM( fse3t_n))) =', z_tmax490 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 492 491 ! 493 492 zht(:,:) = 0.0_wp 494 493 DO jk = 1, jpkm1 495 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)494 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 496 495 END DO 497 496 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 498 497 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 499 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM( fse3t_a))) =', z_tmax498 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 500 499 ! 501 500 zht(:,:) = 0.0_wp 502 501 DO jk = 1, jpkm1 503 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)502 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 504 503 END DO 505 504 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 506 505 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM( fse3t_b))) =', z_tmax506 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 508 507 ! 509 508 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) … … 524 523 ! *********************************** ! 525 524 526 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )527 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )525 CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 526 CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 528 527 529 528 ! *********************************** ! … … 531 530 ! *********************************** ! 532 531 533 hu_a(:,:) = 0._wp ! Ocean depth at U-points534 hv_a(:,:) = 0._wp ! Ocean depth at V-points535 DO jk = 1, jpkm1536 hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)537 hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)532 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 533 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 534 DO jk = 2, jpkm1 535 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 536 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 538 537 END DO 539 538 ! ! Inverse of the local depth 540 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 541 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 542 543 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 544 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 545 539 !!gm BUG ? don't understand the use of umask_i here ..... 540 r1_hu_a(:,:) = umask_i(:,:) / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) 541 r1_hv_a(:,:) = vmask_i(:,:) / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) 542 ! 543 CALL wrk_dealloc( jpi,jpj, zht, z_scale, zwu, zwv, zhdiv ) 544 CALL wrk_dealloc( jpi,jpj,jpk, ze3t ) 545 ! 546 546 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt') 547 547 ! 548 548 END SUBROUTINE dom_vvl_sf_nxt 549 549 … … 561 561 !! - recompute depths and water height fields 562 562 !! 563 !! ** Action : - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step563 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step 564 564 !! - Recompute: 565 !! fse3(u/v)_b566 !! fse3w_n567 !! fse3(u/v)w_b568 !! fse3(u/v)w_n569 !! fsdept_n, fsdepw_n and fsde3w_n565 !! e3(u/v)_b 566 !! e3w_n 567 !! e3(u/v)w_b 568 !! e3(u/v)w_n 569 !! gdept_n, gdepw_n and gde3w_n 570 570 !! h(u/v) and h(u/v)r 571 571 !! … … 573 573 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 574 574 !!---------------------------------------------------------------------- 575 !! * Arguments 576 INTEGER, INTENT( in ) :: kt ! time step 577 !! * Local declarations 578 INTEGER :: ji,jj,jk ! dummy loop indices 579 REAL(wp) :: zcoef 580 !!---------------------------------------------------------------------- 581 575 INTEGER, INTENT( in ) :: kt ! time step 576 ! 577 INTEGER :: ji, jj, jk ! dummy loop indices 578 REAL(wp) :: zcoef ! local scalar 579 !!---------------------------------------------------------------------- 580 ! 581 IF( ln_linssh ) RETURN ! No calculation in linear free surface 582 ! 582 583 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 583 584 ! … … 590 591 ! Time filter and swap of scale factors 591 592 ! ===================================== 592 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.593 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 593 594 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 594 595 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 600 601 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 601 602 ENDIF 602 fsdept_b(:,:,:) = fsdept_n(:,:,:)603 fsdepw_b(:,:,:) = fsdepw_n(:,:,:)604 605 fse3t_n(:,:,:) = fse3t_a(:,:,:)606 fse3u_n(:,:,:) = fse3u_a(:,:,:)607 fse3v_n(:,:,:) = fse3v_a(:,:,:)603 gdept_b(:,:,:) = gdept_n(:,:,:) 604 gdepw_b(:,:,:) = gdepw_n(:,:,:) 605 606 e3t_n(:,:,:) = e3t_a(:,:,:) 607 e3u_n(:,:,:) = e3u_a(:,:,:) 608 e3v_n(:,:,:) = e3v_a(:,:,:) 608 609 609 610 ! Compute all missing vertical scale factor and depths … … 611 612 ! Horizontal scale factor interpolations 612 613 ! -------------------------------------- 613 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt614 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 614 615 ! - JC - hu_b, hv_b, hur_b, hvr_b also 615 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 616 617 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 618 616 619 ! Vertical scale factor interpolations 617 ! ------------------------------------ 618 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 619 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 620 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 621 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 622 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 623 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 624 ! t- and w- points depth 625 ! ---------------------- 626 ! set the isf depth as it is in the initial step 627 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 628 fsdepw_n(:,:,1) = 0.0_wp 629 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 630 620 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 621 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 622 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 623 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' ) 624 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 625 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 626 627 ! t- and w- points depth (set the isf depth as it is in the initial step) 628 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 629 gdepw_n(:,:,1) = 0.0_wp 630 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 631 631 DO jk = 2, jpk 632 632 DO jj = 1,jpj … … 635 635 ! 1 for jk = mikt 636 636 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 637 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)638 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &639 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))640 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)637 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 638 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 639 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 640 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 641 641 END DO 642 642 END DO 643 643 END DO 644 644 645 ! Local depth and Inverse of the local depth of the water column at u- and v- points 646 ! ---------------------------------------------------------------------------------- 647 hu (:,:) = hu_a (:,:) 648 hv (:,:) = hv_a (:,:) 649 650 ! Inverse of the local depth 651 hur(:,:) = hur_a(:,:) 652 hvr(:,:) = hvr_a(:,:) 653 654 ! Local depth of the water column at t- points 655 ! -------------------------------------------- 656 ht(:,:) = 0. 657 DO jk = 1, jpkm1 658 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 645 ! Local depth and Inverse of the local depth of the water 646 ! ------------------------------------------------------- 647 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 648 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 649 ! 650 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 651 DO jk = 2, jpkm1 652 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 659 653 END DO 660 654 661 655 ! Write outputs 662 656 ! ============= 663 CALL iom_put( "e3t" , fse3t_n(:,:,:) )664 CALL iom_put( "e3u" , fse3u_n(:,:,:) )665 CALL iom_put( "e3v" , fse3v_n(:,:,:) )666 CALL iom_put( "e3w" , fse3w_n(:,:,:) )667 CALL iom_put( "tpt_dep" , fsde3w_n(:,:,:) )657 CALL iom_put( "e3t", e3t_n(:,:,:) ) 658 CALL iom_put( "e3u", e3u_n(:,:,:) ) 659 CALL iom_put( "e3v", e3v_n(:,:,:) ) 660 CALL iom_put( "e3w", e3w_n(:,:,:) ) 661 CALL iom_put( "tpt_dep", gde3w_n(:,:,:) ) 668 662 IF( iom_use("e3tdef") ) & 669 CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100* tmask(:,:,:) ) ** 2 )663 CALL iom_put( "e3tdef", ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100. * tmask(:,:,:) ) ** 2 ) 670 664 671 665 ! write restart file 672 666 ! ================== 673 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )674 ! 675 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')676 667 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 668 ! 669 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 670 ! 677 671 END SUBROUTINE dom_vvl_sf_swp 678 672 … … 696 690 !!---------------------------------------------------------------------- 697 691 ! 698 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol')692 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 699 693 ! 700 694 SELECT CASE ( pout ) !== type of interpolation ==! … … 770 764 END SELECT 771 765 ! 772 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol')766 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 773 767 ! 774 768 END SUBROUTINE dom_vvl_interpol … … 801 795 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn ) 802 796 ! 803 id1 = iom_varid( numror, ' fse3t_b', ldstop = .FALSE. )804 id2 = iom_varid( numror, ' fse3t_n', ldstop = .FALSE. )797 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 798 id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 805 799 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 806 800 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) … … 810 804 ! ! --------- ! 811 805 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 812 CALL iom_get( numror, jpdom_autoglo, ' fse3t_b', fse3t_b(:,:,:) )813 CALL iom_get( numror, jpdom_autoglo, ' fse3t_n', fse3t_n(:,:,:) )806 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 807 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 814 808 ! needed to restart if land processor not computed 815 IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'809 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 816 810 WHERE ( tmask(:,:,:) == 0.0_wp ) 817 fse3t_n(:,:,:) = e3t_0(:,:,:)818 fse3t_b(:,:,:) = e3t_0(:,:,:)811 e3t_n(:,:,:) = e3t_0(:,:,:) 812 e3t_b(:,:,:) = e3t_0(:,:,:) 819 813 END WHERE 820 814 IF( neuler == 0 ) THEN 821 fse3t_b(:,:,:) = fse3t_n(:,:,:)815 e3t_b(:,:,:) = e3t_n(:,:,:) 822 816 ENDIF 823 817 ELSE IF( id1 > 0 ) THEN 824 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'825 IF(lwp) write(numout,*) ' fse3t_n set equal to fse3t_b.'818 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 819 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 826 820 IF(lwp) write(numout,*) 'neuler is forced to 0' 827 CALL iom_get( numror, jpdom_autoglo, ' fse3t_b', fse3t_b(:,:,:) )828 fse3t_n(:,:,:) = fse3t_b(:,:,:)821 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 822 e3t_n(:,:,:) = e3t_b(:,:,:) 829 823 neuler = 0 830 824 ELSE IF( id2 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'832 IF(lwp) write(numout,*) ' fse3t_b set equal to fse3t_n.'825 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 826 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 833 827 IF(lwp) write(numout,*) 'neuler is forced to 0' 834 CALL iom_get( numror, jpdom_autoglo, ' fse3t_n', fse3t_n(:,:,:) )835 fse3t_b(:,:,:) = fse3t_n(:,:,:)828 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 829 e3t_b(:,:,:) = e3t_n(:,:,:) 836 830 neuler = 0 837 831 ELSE 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'832 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 839 833 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 840 834 IF(lwp) write(numout,*) 'neuler is forced to 0' 841 DO jk =1,jpk842 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &843 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)&844 & + e3t_0(:,:,jk)* (1._wp -tmask(:,:,jk))835 DO jk = 1, jpk 836 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 837 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 838 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 845 839 END DO 846 fse3t_b(:,:,:) = fse3t_n(:,:,:)840 e3t_b(:,:,:) = e3t_n(:,:,:) 847 841 neuler = 0 848 842 ENDIF … … 875 869 ! 876 870 ELSE !* Initialize at "rest" 877 fse3t_b(:,:,:) = e3t_0(:,:,:)878 fse3t_n(:,:,:) = e3t_0(:,:,:)871 e3t_b(:,:,:) = e3t_0(:,:,:) 872 e3t_n(:,:,:) = e3t_0(:,:,:) 879 873 sshn(:,:) = 0.0_wp 880 874 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN … … 891 885 ! ! all cases ! 892 886 ! ! --------- ! 893 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_b', fse3t_b(:,:,:) )894 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_n', fse3t_n(:,:,:) )887 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) ) 888 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) ) 895 889 ! ! ----------------------- ! 896 890 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 975 969 ! 976 970 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 977 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0)CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )971 IF( .NOT. ln_vvl_zstar .AND. nn_isf /= 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 978 972 ! 979 973 IF(lwp) THEN ! Print the choice … … 989 983 ! 990 984 #if defined key_agrif 991 IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )985 IF(.NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' ) 992 986 #endif 993 987 ! … … 996 990 !!====================================================================== 997 991 END MODULE domvvl 998 999 1000 -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5836 r6060 2 2 !!============================================================================== 3 3 !! *** MODULE domzgr *** 4 !! Ocean initialization : domain initialization4 !! Ocean domain : definition of the vertical coordinate system 5 5 !!============================================================================== 6 6 !! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate … … 38 38 USE closea ! closed seas 39 39 USE c1d ! 1D vertical configuration 40 ! 40 41 USE in_out_manager ! I/O manager 41 42 USE iom ! I/O library … … 73 74 74 75 !! * Substitutions 75 # include "domzgr_substitute.h90"76 76 # include "vectopt_loop_substitute.h90" 77 77 !!---------------------------------------------------------------------- … … 102 102 INTEGER :: ios 103 103 ! 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 105 105 !!---------------------------------------------------------------------- 106 106 ! … … 120 120 WRITE(numout,*) 'dom_zgr : vertical coordinate' 121 121 WRITE(numout,*) '~~~~~~~' 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 127 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 127 128 ENDIF 129 130 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 128 131 129 132 ioptio = 0 ! Check Vertical coordinate options … … 155 158 ! 156 159 IF( nprint == 1 .AND. lwp ) THEN 157 WRITE(numout,*) ' MIN val mbathy ', MINVAL(mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )160 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 158 161 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 159 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )160 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL(e3f_0(:,:,:) ), &161 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL(e3v_0(:,:,:) ), &162 & ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL(e3vw_0(:,:,:)), &163 & ' w ', MINVAL(e3w_0(:,:,:) )162 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) 163 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), & 164 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), & 165 & ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), & 166 & ' w ', MINVAL( e3w_0(:,:,:) ) 164 167 165 168 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 166 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )167 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL(e3f_0(:,:,:) ), &168 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL(e3v_0(:,:,:) ), &169 & ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),&170 & ' w ', MAXVAL(e3w_0(:,:,:) )169 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) ) 170 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), & 171 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), & 172 & ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), & 173 & ' w ', MAXVAL( e3w_0(:,:,:) ) 171 174 ENDIF 172 175 ! … … 674 677 !! - update bathy : meter bathymetry (in meters) 675 678 !!---------------------------------------------------------------------- 676 !!677 679 INTEGER :: ji, jj, jl ! dummy loop indices 678 680 INTEGER :: icompt, ibtest, ikmax ! temporary integers 679 681 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 680 681 682 !!---------------------------------------------------------------------- 682 683 ! … … 775 776 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 776 777 ENDIF 777 778 IF( lwp .AND. nprint == 1 ) THEN ! control print779 WRITE(numout,*)780 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '781 WRITE(numout,*) ' ------------------'782 CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )783 WRITE(numout,*)784 ENDIF785 778 ! 786 779 CALL wrk_dealloc( jpi, jpj, zbathy ) … … 803 796 !! (min value = 1 over land) 804 797 !!---------------------------------------------------------------------- 805 !!806 798 INTEGER :: ji, jj ! dummy loop indices 807 799 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk … … 835 827 END SUBROUTINE zgr_bot_level 836 828 837 SUBROUTINE zgr_top_level 829 830 SUBROUTINE zgr_top_level 838 831 !!---------------------------------------------------------------------- 839 832 !! *** ROUTINE zgr_bot_level *** … … 847 840 !! (min value = 1) 848 841 !!---------------------------------------------------------------------- 849 !!850 842 INTEGER :: ji, jj ! dummy loop indices 851 843 REAL(wp), POINTER, DIMENSION(:,:) :: zmik … … 881 873 END SUBROUTINE zgr_top_level 882 874 875 883 876 SUBROUTINE zgr_zco 884 877 !!---------------------------------------------------------------------- 885 878 !! *** ROUTINE zgr_zco *** 886 879 !! 887 !! ** Purpose : define the z-coordinate system880 !! ** Purpose : define the reference z-coordinate system 888 881 !! 889 882 !! ** Method : set 3D coord. arrays to reference 1D array … … 895 888 ! 896 889 DO jk = 1, jpk 897 gdept_0 898 gdepw_0 899 gde p3w_0(:,:,jk) = gdepw_1d(jk)900 e3t_0 901 e3u_0 902 e3v_0 903 e3f_0 904 e3w_0 905 e3uw_0 906 e3vw_0 890 gdept_0(:,:,jk) = gdept_1d(jk) 891 gdepw_0(:,:,jk) = gdepw_1d(jk) 892 gde3w_0(:,:,jk) = gdepw_1d(jk) 893 e3t_0 (:,:,jk) = e3t_1d (jk) 894 e3u_0 (:,:,jk) = e3t_1d (jk) 895 e3v_0 (:,:,jk) = e3t_1d (jk) 896 e3f_0 (:,:,jk) = e3t_1d (jk) 897 e3w_0 (:,:,jk) = e3w_1d (jk) 898 e3uw_0 (:,:,jk) = e3w_1d (jk) 899 e3vw_0 (:,:,jk) = e3w_1d (jk) 907 900 END DO 908 901 ! … … 917 910 !! 918 911 !! ** Purpose : the depth and vertical scale factor in partial step 919 !! z-coordinate case912 !! reference z-coordinate case 920 913 !! 921 914 !! ** Method : Partial steps : computes the 3D vertical scale factors … … 957 950 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 958 951 !!---------------------------------------------------------------------- 959 !!960 952 INTEGER :: ji, jj, jk ! dummy loop indices 961 953 INTEGER :: ik, it, ikb, ikt ! temporary integers 962 LOGICAL :: ll_print ! Allow control print for debugging963 954 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 964 955 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t … … 971 962 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 972 963 ! 973 CALL wrk_alloc( jpi, jpj, jpk,zprt )964 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 974 965 ! 975 966 IF(lwp) WRITE(numout,*) … … 977 968 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 978 969 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 979 980 ll_print = .FALSE. ! Local variable for debugging981 982 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth983 WRITE(numout,*)984 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'985 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )986 ENDIF987 988 970 989 971 ! bathymetry in level (from bathy_meter) … … 1196 1178 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1197 1179 1198 ! Compute gde p3w_0 (vertical sum of e3w)1180 ! Compute gde3w_0 (vertical sum of e3w) 1199 1181 IF ( ln_isfcav ) THEN ! if cavity 1200 WHERE (misfdep == 0)misfdep = 11182 WHERE( misfdep == 0 ) misfdep = 1 1201 1183 DO jj = 1,jpj 1202 1184 DO ji = 1,jpi 1203 gde p3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1185 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1204 1186 DO jk = 2, misfdep(ji,jj) 1205 gde p3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1187 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1206 1188 END DO 1207 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))1189 IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1208 1190 DO jk = misfdep(ji,jj) + 1, jpk 1209 gde p3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1191 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1210 1192 END DO 1211 1193 END DO 1212 1194 END DO 1213 1195 ELSE ! no cavity 1214 gde p3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1196 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1215 1197 DO jk = 2, jpk 1216 gde p3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)1198 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1217 1199 END DO 1218 1200 END IF 1219 ! ! ================= ! 1220 IF(lwp .AND. ll_print) THEN ! Control print ! 1221 ! ! ================= ! 1222 DO jj = 1,jpj 1223 DO ji = 1, jpi 1224 ik = MAX( mbathy(ji,jj), 1 ) 1225 zprt(ji,jj,1) = e3t_0 (ji,jj,ik) 1226 zprt(ji,jj,2) = e3w_0 (ji,jj,ik) 1227 zprt(ji,jj,3) = e3u_0 (ji,jj,ik) 1228 zprt(ji,jj,4) = e3v_0 (ji,jj,ik) 1229 zprt(ji,jj,5) = e3f_0 (ji,jj,ik) 1230 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 1231 END DO 1232 END DO 1233 WRITE(numout,*) 1234 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1235 WRITE(numout,*) 1236 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1237 WRITE(numout,*) 1238 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1239 WRITE(numout,*) 1240 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1241 WRITE(numout,*) 1242 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1243 WRITE(numout,*) 1244 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1245 ENDIF 1246 ! 1247 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1201 ! 1202 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1248 1203 ! 1249 1204 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1250 1205 ! 1251 1206 END SUBROUTINE zgr_zps 1207 1252 1208 1253 1209 SUBROUTINE zgr_isf … … 1265 1221 !! - bathy and isfdraft are modified 1266 1222 !!---------------------------------------------------------------------- 1267 !!1268 1223 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1269 1224 INTEGER :: ik, it ! temporary integers 1270 1225 INTEGER :: id, jd, nprocd 1271 1226 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 1272 LOGICAL :: ll_print ! Allow control print for debugging1273 1227 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1274 1228 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t … … 1281 1235 !!--------------------------------------------------------------------- 1282 1236 ! 1283 IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1284 ! 1285 CALL wrk_alloc( jpi, jpj,zbathy, zmask, zrisfdep)1286 CALL wrk_alloc( jpi, jpj,zmisfdep, zmbathy )1237 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1238 ! 1239 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1240 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1287 1241 1288 1242 … … 1300 1254 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1301 1255 END DO 1302 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT.0._wp)1303 risfdep(:,:) = 0. ;misfdep(:,:) = 11256 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) > 0._wp) 1257 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1304 1258 END WHERE 1305 1259 … … 1308 1262 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1309 1263 DO jl = 1, 10 1310 WHERE (bathy(:,:) .EQ.risfdep(:,:) )1264 WHERE (bathy(:,:) == risfdep(:,:) ) 1311 1265 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1312 1266 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp … … 1315 1269 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1316 1270 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1317 END WHERE1271 END WHERE 1318 1272 IF( lk_mpp ) THEN 1319 1273 zbathy(:,:) = FLOAT( misfdep(:,:) ) … … 1360 1314 ! find the minimum change option: 1361 1315 ! test bathy 1362 IF (risfdep(ji,jj) .GT.1) THEN1316 IF (risfdep(ji,jj) > 1) THEN 1363 1317 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1364 1318 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) … … 1366 1320 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1367 1321 1368 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT.misfdep(ji,jj)) THEN1322 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1369 1323 IF (zbathydiff .LE. zrisfdepdiff) THEN 1370 1324 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) … … 1752 1706 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1753 1707 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1754 1755 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf')1756 1708 ! 1709 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1710 ! 1757 1711 END SUBROUTINE 1712 1758 1713 1759 1714 SUBROUTINE zgr_sco … … 1801 1756 !! 1802 1757 !!---------------------------------------------------------------------- 1803 !1804 1758 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1805 1759 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers … … 1810 1764 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1811 1765 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1812 1766 !! 1813 1767 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1814 1768 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1815 1769 !!---------------------------------------------------------------------- 1816 1770 ! 1817 1771 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1818 1772 ! 1819 CALL wrk_alloc( jpi, jpj,zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )1773 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1820 1774 ! 1821 1775 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 1876 1830 DO jj = 1, jpj 1877 1831 DO ji = 1, jpi 1878 IF( bathy(ji,jj) == 0._wp ) THEN 1879 iip1 = MIN( ji+1, jpi ) 1880 ijp1 = MIN( jj+1, jpj ) 1881 iim1 = MAX( ji-1, 1 ) 1882 ijm1 = MAX( jj-1, 1 ) 1883 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1884 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1885 zenv(ji,jj) = rn_sbot_min 1886 ENDIF 1832 IF( bathy(ji,jj) == 0._wp ) THEN 1833 iip1 = MIN( ji+1, jpi ) 1834 ijp1 = MIN( jj+1, jpj ) 1835 iim1 = MAX( ji-1, 1 ) 1836 ijm1 = MAX( jj-1, 1 ) 1837 !!gm BUG fix see ticket #1617 1838 IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) & 1839 & + bathy(iim1,jj ) + bathy(iip1,jj ) & 1840 & + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) zenv(ji,jj) = rn_sbot_min 1841 !!gm 1842 !!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + & 1843 !!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1844 !!gm zenv(ji,jj) = rn_sbot_min 1845 !!gm ENDIF 1846 !!gm end 1887 1847 ENDIF 1888 1848 END DO … … 1975 1935 ENDIF 1976 1936 ! 1977 IF(lwp) THEN ! Control print1978 WRITE(numout,*)1979 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'1980 WRITE(numout,*)1981 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )1982 IF( nprint == 1 ) THEN1983 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )1984 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )1985 ENDIF1986 ENDIF1987 1988 1937 ! ! ============================== 1989 1938 ! ! hbatu, hbatv, hbatf fields … … 2080 2029 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 2081 2030 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2082 2083 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2084 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2085 ! 2086 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2087 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2088 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2089 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2090 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2091 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2092 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2031 ! 2032 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2033 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2034 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2035 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2036 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2037 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2038 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2093 2039 2094 2040 #if defined key_agrif 2095 ! Ensure meaningful vertical scale factors in ghost lines/columns 2096 IF( .NOT. Agrif_Root() ) THEN 2097 ! 2098 IF((nbondi == -1).OR.(nbondi == 2)) THEN 2099 e3u_0(1,:,:) = e3u_0(2,:,:) 2100 ENDIF 2101 ! 2102 IF((nbondi == 1).OR.(nbondi == 2)) THEN 2103 e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) 2104 ENDIF 2105 ! 2106 IF((nbondj == -1).OR.(nbondj == 2)) THEN 2107 e3v_0(:,1,:) = e3v_0(:,2,:) 2108 ENDIF 2109 ! 2110 IF((nbondj == 1).OR.(nbondj == 2)) THEN 2111 e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) 2112 ENDIF 2113 ! 2114 ENDIF 2041 IF( .NOT. Agrif_Root() ) THEN ! Ensure meaningful vertical scale factors in ghost lines/columns 2042 IF( nbondi == -1 .OR. nbondi == 2 ) e3u_0( 1 , : ,:) = e3u_0( 2 , : ,:) 2043 IF( nbondi == 1 .OR. nbondi == 2 ) e3u_0(nlci-1, : ,:) = e3u_0(nlci-2, : ,:) 2044 IF( nbondj == -1 .OR. nbondj == 2 ) e3v_0( : , 1 ,:) = e3v_0( : , 2 ,:) 2045 IF( nbondj == 1 .OR. nbondj == 2 ) e3v_0( : ,nlcj-1,:) = e3v_0( : ,nlcj-2,:) 2046 ENDIF 2115 2047 #endif 2116 2048 2117 fsdept(:,:,:) = gdept_0 (:,:,:) 2118 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2119 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2120 fse3t (:,:,:) = e3t_0 (:,:,:) 2121 fse3u (:,:,:) = e3u_0 (:,:,:) 2122 fse3v (:,:,:) = e3v_0 (:,:,:) 2123 fse3f (:,:,:) = e3f_0 (:,:,:) 2124 fse3w (:,:,:) = e3w_0 (:,:,:) 2125 fse3uw(:,:,:) = e3uw_0 (:,:,:) 2126 fse3vw(:,:,:) = e3vw_0 (:,:,:) 2049 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 2050 !!gm and only that !!!!! 2051 !!gm THIS should be removed from here ! 2052 gdept_n(:,:,:) = gdept_0(:,:,:) 2053 gdepw_n(:,:,:) = gdepw_0(:,:,:) 2054 gde3w_n(:,:,:) = gde3w_0(:,:,:) 2055 e3t_n (:,:,:) = e3t_0 (:,:,:) 2056 e3u_n (:,:,:) = e3u_0 (:,:,:) 2057 e3v_n (:,:,:) = e3v_0 (:,:,:) 2058 e3f_n (:,:,:) = e3f_0 (:,:,:) 2059 e3w_n (:,:,:) = e3w_0 (:,:,:) 2060 e3uw_n (:,:,:) = e3uw_0 (:,:,:) 2061 e3vw_n (:,:,:) = e3vw_0 (:,:,:) 2062 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 2063 !! gm end 2127 2064 !! 2128 2065 ! HYBRID : … … 2130 2067 DO ji = 1, jpi 2131 2068 DO jk = 1, jpkm1 2132 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )2133 END DO 2134 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 02069 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2070 END DO 2071 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2135 2072 END DO 2136 2073 END DO … … 2141 2078 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2142 2079 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2143 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde p3w_0(:,:,:) )2144 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 2145 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 2146 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 2080 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 2081 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2082 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 2083 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 2147 2084 & ' w ', MINVAL( e3w_0 (:,:,:) ) 2148 2085 2149 2086 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2150 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde p3w_0(:,:,:) )2151 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 2152 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 2153 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 2087 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 2088 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2089 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 2090 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 2154 2091 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 2155 2092 ENDIF … … 2183 2120 END DO 2184 2121 ENDIF 2185 2186 !================================================================================2187 ! check the coordinate makes sense2188 !================================================================================2122 ! 2123 !================================================================================ 2124 ! check the coordinate makes sense 2125 !================================================================================ 2189 2126 DO ji = 1, jpi 2190 2127 DO jj = 1, jpj 2191 2128 ! 2192 2129 IF( hbatt(ji,jj) > 0._wp) THEN 2193 2130 DO jk = 1, mbathy(ji,jj) 2194 2131 ! check coordinate is monotonically increasing 2195 IF ( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN2132 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 2196 2133 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2197 2134 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2198 WRITE(numout,*) 'e3w', fse3w(ji,jj,:)2199 WRITE(numout,*) 'e3t', fse3t(ji,jj,:)2135 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2136 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 2200 2137 CALL ctl_stop( ctmp1 ) 2201 2138 ENDIF 2202 2139 ! and check it has never gone negative 2203 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN2140 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 2204 2141 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2205 2142 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2206 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2207 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2143 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2144 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2208 2145 CALL ctl_stop( ctmp1 ) 2209 2146 ENDIF 2210 2147 ! and check it never exceeds the total depth 2211 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN2148 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2212 2149 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2213 2150 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2214 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2151 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2215 2152 CALL ctl_stop( ctmp1 ) 2216 2153 ENDIF 2217 2154 END DO 2218 2155 ! 2219 2156 DO jk = 1, mbathy(ji,jj)-1 2220 2157 ! and check it never exceeds the total depth 2221 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN2158 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2222 2159 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2223 2160 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2224 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2161 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2225 2162 CALL ctl_stop( ctmp1 ) 2226 2163 ENDIF 2227 2164 END DO 2228 2229 2165 ENDIF 2230 2231 2166 END DO 2232 2167 END DO … … 2238 2173 END SUBROUTINE zgr_sco 2239 2174 2240 !!====================================================================== 2175 2241 2176 SUBROUTINE s_sh94() 2242 2243 2177 !!---------------------------------------------------------------------- 2244 2178 !! *** ROUTINE s_sh94 *** … … 2251 2185 !! Reference : Song and Haidvogel 1994. 2252 2186 !!---------------------------------------------------------------------- 2253 !2254 2187 INTEGER :: ji, jj, jk ! dummy loop argument 2255 2188 REAL(wp) :: zcoeft, zcoefw ! temporary scalars … … 2257 2190 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2258 2191 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2259 2260 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2261 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2192 !!---------------------------------------------------------------------- 2193 2194 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2195 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2262 2196 2263 2197 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp … … 2265 2199 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2266 2200 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2267 2201 ! 2268 2202 DO ji = 1, jpi 2269 2203 DO jj = 1, jpj 2270 2204 ! 2271 2205 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 2272 2206 DO jk = 1, jpk … … 2297 2231 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2298 2232 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2299 gdept_0 2300 gdepw_0 2301 gde p3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )2233 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2234 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2235 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 2302 2236 END DO 2303 2237 ! … … 2331 2265 END DO 2332 2266 END DO 2333 2334 CALL wrk_dealloc( jpi, jpj, jpk,z_gsigw3, z_gsigt3, z_gsi3w3 )2335 CALL wrk_dealloc( jpi, jpj, jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )2336 2267 ! 2268 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2269 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2270 ! 2337 2271 END SUBROUTINE s_sh94 2338 2272 2273 2339 2274 SUBROUTINE s_sf12 2340 2341 2275 !!---------------------------------------------------------------------- 2342 2276 !! *** ROUTINE s_sf12 *** … … 2354 2288 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 2355 2289 !!---------------------------------------------------------------------- 2356 !2357 2290 INTEGER :: ji, jj, jk ! dummy loop argument 2358 2291 REAL(wp) :: zsmth ! smoothing around critical depth … … 2361 2294 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2362 2295 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2363 2296 !!---------------------------------------------------------------------- 2364 2297 ! 2365 2298 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) … … 2425 2358 2426 2359 DO jk = 1, jpk 2427 gdept_0 2428 gdepw_0 2429 gde p3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)2360 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2361 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2362 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 2430 2363 END DO 2431 2364 … … 2454 2387 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 2455 2388 ! 2456 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)2389 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 2457 2390 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 2458 2391 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) … … 2467 2400 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 2468 2401 ! 2469 ! ! ============= 2470 2471 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2472 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2473 2402 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2403 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2404 ! 2474 2405 END SUBROUTINE s_sf12 2475 2406 2407 2476 2408 SUBROUTINE s_tanh() 2477 2478 2409 !!---------------------------------------------------------------------- 2479 2410 !! *** ROUTINE s_tanh*** … … 2485 2416 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 2486 2417 !!---------------------------------------------------------------------- 2487 2488 INTEGER :: ji, jj, jk ! dummy loop argument 2418 INTEGER :: ji, jj, jk ! dummy loop argument 2489 2419 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2490 2491 2420 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 2492 2421 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 2493 2494 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2495 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2422 !!---------------------------------------------------------------------- 2423 2424 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2425 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2496 2426 2497 2427 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp … … 2523 2453 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2524 2454 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2525 gdept_0 2526 gdepw_0 2527 gde p3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )2455 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2456 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2457 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2528 2458 END DO 2529 2459 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 2542 2472 END DO 2543 2473 END DO 2544 2545 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w)2546 CALL wrk_dealloc( jpk, z_esigt, z_esigw)2547 2474 ! 2475 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2476 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2477 ! 2548 2478 END SUBROUTINE s_tanh 2479 2549 2480 2550 2481 FUNCTION fssig( pk ) RESULT( pf ) … … 2618 2549 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 2619 2550 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 2620 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth2621 REAL(wp), INTENT(in ) :: pzs ! surface box depth2622 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter2623 REAL(wp) :: za1,za2,za3 ! local variables2624 REAL(wp) :: zn1,zn2 ! local variables2625 REAL(wp) :: za,zb,zx ! local variables2626 integer :: jk2627 !!----------------------------------------------------------------------2628 ! 2629 2630 zn1 = 1. /(jpk-1.)2631 zn2 = 1. - zn12632 2551 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 2552 REAL(wp), INTENT(in ) :: pzs ! surface box depth 2553 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 2554 ! 2555 INTEGER :: jk ! dummy loop index 2556 REAL(wp) :: za1,za2,za3 ! local scalar 2557 REAL(wp) :: zn1,zn2 ! - - 2558 REAL(wp) :: za,zb,zx ! - - 2559 !!---------------------------------------------------------------------- 2560 ! 2561 zn1 = 1._wp / REAL( jpkm1, wp ) 2562 zn2 = 1._wp - zn1 2563 ! 2633 2564 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 2634 2565 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 2635 2566 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2636 2567 ! 2637 2568 za = pzb - za3*(pzs-za1)-za2 2638 2569 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 2639 2570 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 2640 2571 zx = 1.0_wp-za/2.0_wp-zb 2641 2572 ! 2642 2573 DO jk = 1, jpk 2643 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + &2644 &zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &2645 &(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )2574 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 2575 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 2576 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 2646 2577 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 2647 ENDDO 2648 2578 END DO 2649 2579 ! 2650 2580 END FUNCTION fgamma -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r5836 r6060 35 35 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsd ! structure of input SST (file informations, fields read) 36 36 37 !! * Substitutions38 # include "domzgr_substitute.h90"39 37 !!---------------------------------------------------------------------- 40 38 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 250 248 ENDIF 251 249 ! 252 IF( lwp .AND. kt == nit000 ) THEN253 WRITE(numout,*) ' temperature Levitus '254 WRITE(numout,*)255 WRITE(numout,*)' level = 1'256 CALL prihre( ptsd(:,:,1 ,jp_tem), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )257 WRITE(numout,*)' level = ', jpk/2258 CALL prihre( ptsd(:,:,jpk/2,jp_tem), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )259 WRITE(numout,*)' level = ', jpkm1260 CALL prihre( ptsd(:,:,jpkm1,jp_tem), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )261 WRITE(numout,*)262 WRITE(numout,*) ' salinity Levitus '263 WRITE(numout,*)264 WRITE(numout,*)' level = 1'265 CALL prihre( ptsd(:,:,1 ,jp_sal), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )266 WRITE(numout,*)' level = ', jpk/2267 CALL prihre( ptsd(:,:,jpk/2,jp_sal), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )268 WRITE(numout,*)' level = ', jpkm1269 CALL prihre( ptsd(:,:,jpkm1,jp_sal), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )270 WRITE(numout,*)271 ENDIF272 !273 250 IF( .NOT.ln_tsd_tradmp ) THEN !== deallocate T & S structure ==! 274 251 ! (data used only for initialisation) -
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5930 r6060 49 49 50 50 !! * Substitutions 51 # include "domzgr_substitute.h90"52 51 # include "vectopt_loop_substitute.h90" 53 52 !!---------------------------------------------------------------------- … … 120 119 ENDIF 121 120 ENDIF 122 ! 123 ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 124 IF( lk_vvl ) THEN 121 ! 122 !!gm This is to be changed !!!! 123 ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here 124 IF( .NOT.ln_linssh ) THEN 125 125 DO jk = 1, jpk 126 fse3t_b(:,:,jk) = fse3t_n(:,:,jk)126 e3t_b(:,:,jk) = e3t_n(:,:,jk) 127 127 END DO 128 128 ENDIF 129 !!gm 129 130 ! 130 131 ENDIF 131 !132 132 ! 133 133 ! Initialize "now" and "before" barotropic velocities: 134 ! Do it whatever the free surface method, these arrays 135 ! being eventually used 136 ! 134 ! Do it whatever the free surface method, these arrays being eventually used 137 135 ! 138 136 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 139 137 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 140 138 ! 139 !!gm the use of umsak & vmask is not necessary belox as un, vn, ub, vb are always masked 141 140 DO jk = 1, jpkm1 142 141 DO jj = 1, jpj 143 142 DO ji = 1, jpi 144 un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)145 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)143 un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 144 vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 146 145 ! 147 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)148 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)146 ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 147 vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 149 148 END DO 150 149 END DO 151 150 END DO 152 151 ! 153 un_b(:,:) = un_b(:,:) * hur (:,:) 154 vn_b(:,:) = vn_b(:,:) * hvr (:,:) 155 ! 156 ub_b(:,:) = ub_b(:,:) * hur_b(:,:) 157 vb_b(:,:) = vb_b(:,:) * hvr_b(:,:) 158 ! 152 un_b(:,:) = un_b(:,:) * r1_hu_n(:,:) 153 vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:) 154 ! 155 ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 156 vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 159 157 ! 160 158 IF( nn_timing == 1 ) CALL timing_stop('istate_init') … … 184 182 ! 185 183 DO jk = 1, jpk 186 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH(( fsdept(:,:,jk)-80.)/30.) ) &187 & + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.) ) * tmask(:,:,jk)184 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((gdept_n(:,:,jk)-80.)/30.) ) & 185 & + 10. * ( 5000. - gdept_n(:,:,jk) ) /5000.) ) * tmask(:,:,jk) 188 186 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 189 187 END DO … … 238 236 ! 239 237 DO jk = 1, jpk 240 tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)238 tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - gdept_n(:,:,jk) / 1000 ) ) * tmask(:,:,jk) 241 239 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 242 240 END DO 243 !244 IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &245 & 1 , jpi , 5 , 1 , jpk , &246 & 1 , 1. , numout )247 241 ! 248 242 ! set salinity field to a constant value … … 314 308 tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb 315 309 ! 316 IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &317 & 1 , jpi , 5 , 1 , jpk , &318 & 1 , 1. , numout )319 !320 310 ! set salinity field to a constant value 321 311 ! -------------------------------------- … … 363 353 DO jj = 1, jpj 364 354 DO ji = 1, jpi 365 tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( ( fsdept(ji,jj,jk) - 400) / 700 ) ) &366 & * (-TANH( (500- fsdept(ji,jj,jk)) / 150 ) + 1) / 2 &367 & + ( 15. * ( 1. - TANH( ( fsdept(ji,jj,jk)-50.) / 1500.) ) &368 & - 1.4 * TANH(( fsdept(ji,jj,jk)-100.) / 100.) &369 & + 7. * (1500. - fsdept(ji,jj,jk)) / 1500. ) &370 & * (-TANH( ( fsdept(ji,jj,jk) - 500) / 150) + 1) / 2355 tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 ) ) & 356 & * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2 & 357 & + ( 15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) ) & 358 & - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.) & 359 & + 7. * (1500. - gdept_n(ji,jj,jk)) / 1500. ) & 360 & * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 371 361 tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 372 362 tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) 373 363 374 tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( ( fsdept(ji,jj,jk) - 305) / 460 ) ) &375 & * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2 &376 & + ( 35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000. &377 & - 1.62 * TANH( ( fsdept(ji,jj,jk) - 60. ) / 650. ) &378 & + 0.2 * TANH( ( fsdept(ji,jj,jk) - 35. ) / 100. ) &379 & + 0.2 * TANH( ( fsdept(ji,jj,jk) - 1000.) / 5000.) ) &380 & * (-TANH(( fsdept(ji,jj,jk) - 500) / 150) + 1) / 2364 tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 ) ) & 365 & * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2 & 366 & + ( 35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000. & 367 & - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60. ) / 650. ) & 368 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 35. ) / 100. ) & 369 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.) ) & 370 & * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 381 371 tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 382 372 tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) … … 451 441 zalfg = 0.5 * grav * rau0 452 442 453 zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value443 zprn(:,:,1) = zalfg * e3w_n(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value 454 444 455 445 DO jk = 2, jpkm1 ! Vertical integration from the surface 456 446 zprn(:,:,jk) = zprn(:,:,jk-1) & 457 & + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )447 & + zalfg * e3w_n(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) ) 458 448 END DO 459 449
Note: See TracChangeset
for help on using the changeset viewer.
