Changeset 6141 for branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM
- Timestamp:
- 2015-12-21T12:38:26+01:00 (8 years ago)
- Location:
- branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM
- Files:
-
- 1 deleted
- 12 edited
- 3 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
r5563 r6141 11 11 !! ! 2004-01 (A.M. Treguier) new calculation based on adatrj 12 12 !! ! 2006-08 (G. Madec) surface module major update 13 !! ! 2015-11 (D. Lea) Allow non-zero initial time of day 13 14 !!---------------------------------------------------------------------- 14 15 … … 20 21 !! 21 22 !! we suppose that the time step is deviding the number of second of in a day 22 !! ---> MOD( rday, rdt tra(1)) == 023 !! ---> MOD( rday, rdt ) == 0 23 24 !! 24 25 !! ----------- WARNING ----------- … … 26 27 !! 27 28 !!---------------------------------------------------------------------- 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 29 USE dom_oce ! ocean space and time domain 30 USE phycst ! physical constants 31 USE in_out_manager ! I/O manager 32 USE iom ! 33 USE ioipsl , ONLY : ymds2ju ! for calendar 34 USE prtctl ! Print control 35 USE trc_oce , ONLY : lk_offline ! offline flag 36 USE timing ! Timing 37 USE restart ! restart 37 38 38 39 IMPLICIT NONE … … 43 44 PUBLIC day_mth ! Needed by TAM 44 45 45 INTEGER, PUBLIC :: nsecd, nsecd05, ndt, ndt05 !(PUBLIC for TAM)46 INTEGER, PUBLIC :: nsecd, nsecd05, ndt, ndt05 !: (PUBLIC for TAM) 46 47 47 48 !!---------------------------------------------------------------------- … … 78 79 & 'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' ) 79 80 ENDIF 80 ! all calendar staff is based on the fact that MOD( rday, rdt tra(1)) == 081 IF( MOD( rday , rdt tra(1)) /= 0. ) CALL ctl_stop( 'the time step must devide the number of second of in a day' )82 IF( MOD( rday , 2. 83 IF( MOD( rdt tra(1), 2.) /= 0. ) CALL ctl_stop( 'the time step (in second) must be an even number' )84 nsecd = NINT(rday 85 nsecd05 = NINT(0.5 * rday 86 ndt = NINT( rdt tra(1))87 ndt05 = NINT(0.5 * rdt tra(1))81 ! all calendar staff is based on the fact that MOD( rday, rdt ) == 0 82 IF( MOD( rday , rdt ) /= 0. ) CALL ctl_stop( 'the time step must devide the number of second of in a day' ) 83 IF( MOD( rday , 2. ) /= 0. ) CALL ctl_stop( 'the number of second of in a day must be an even number' ) 84 IF( MOD( rdt , 2. ) /= 0. ) CALL ctl_stop( 'the time step (in second) must be an even number' ) 85 nsecd = NINT(rday ) 86 nsecd05 = NINT(0.5 * rday ) 87 ndt = NINT( rdt ) 88 ndt05 = NINT(0.5 * rdt ) 88 89 89 90 IF( .NOT. lk_offline ) CALL day_rst( nit000, 'READ' ) … … 95 96 nday = ndastp - (nyear * 10000) - ( nmonth * 100 ) 96 97 97 CALL ymds2ju( nyear, nmonth, nday, 0.0, fjulday ) ! we assume that we start run at 00:00 98 nhour = nn_time0 / 100 99 nminute = ( nn_time0 - nhour * 100 ) 100 101 CALL ymds2ju( nyear, nmonth, nday, nhour*3600._wp+nminute*60._wp, fjulday ) 98 102 IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < 0.1 / rday ) fjulday = REAL(NINT(fjulday),wp) ! avoid truncation error 99 fjulday = fjulday + 1.! move back to the day at nit000 (and not at nit000 - 1)103 IF( nn_time0*3600 - ndt05 .lt. 0 ) fjulday = fjulday + 1. ! move back to the day at nit000 (and not at nit000 - 1) 100 104 101 105 nsec1jan000 = 0 … … 118 122 !compute number of days between last monday and today 119 123 CALL ymds2ju( 1900, 01, 01, 0.0, zjul ) ! compute julian day value of 01.01.1900 (our reference that was a Monday) 120 inbday = NINT(fjulday - zjul) ! compute nb day between 01.01.1900 andcurrent day124 inbday = FLOOR(fjulday - zjul) ! compute nb day between 01.01.1900 and start of current day 121 125 idweek = MOD(inbday, 7) ! compute nb day between last monday and current day 126 IF (idweek .lt. 0) idweek=idweek+7 ! Avoid negative values for dates before 01.01.1900 122 127 123 128 ! number of seconds since the beginning of current year/month/week/day at the middle of the time-step 124 nsec_year = nday_year * nsecd - ndt05 ! 1 time step before the middle of the first time step 125 nsec_month = nday * nsecd - ndt05 ! because day will be called at the beginning of step 126 nsec_week = idweek * nsecd - ndt05 127 nsec_day = nsecd - ndt05 129 IF (nhour*3600+nminute*60-ndt05 .gt. 0) THEN 130 ! 1 timestep before current middle of first time step is still the same day 131 nsec_year = (nday_year-1) * nsecd + nhour*3600+nminute*60 - ndt05 132 nsec_month = (nday-1) * nsecd + nhour*3600+nminute*60 - ndt05 133 ELSE 134 ! 1 time step before the middle of the first time step is the previous day 135 nsec_year = nday_year * nsecd + nhour*3600+nminute*60 - ndt05 136 nsec_month = nday * nsecd + nhour*3600+nminute*60 - ndt05 137 ENDIF 138 nsec_week = idweek * nsecd + nhour*3600+nminute*60 - ndt05 139 nsec_day = nhour*3600+nminute*60 - ndt05 140 IF( nsec_day .lt. 0 ) nsec_day = nsec_day + nsecd 141 IF( nsec_week .lt. 0 ) nsec_week = nsec_week + nsecd*7 128 142 129 143 ! control print 130 IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i8,a,i8)')' =======>> 1/2 time step before the start of the run DATE Y/M/D = ', & 131 & nyear, '/', nmonth, '/', nday, ' nsec_day:', nsec_day, ' nsec_week:', nsec_week 144 IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i8,a,i8,a,i8,a,i8)')' =======>> 1/2 time step before the start of the run DATE Y/M/D = ', & 145 & nyear, '/', nmonth, '/', nday, ' nsec_day:', nsec_day, ' nsec_week:', nsec_week, ' & 146 & nsec_month:', nsec_month , ' nsec_year:' , nsec_year 132 147 133 148 ! Up to now, calendar parameters are related to the end of previous run (nit000-1) … … 223 238 nsec_week = nsec_week + ndt 224 239 nsec_day = nsec_day + ndt 225 adatrj = adatrj + rdt tra(1)/ rday226 fjulday = fjulday + rdt tra(1)/ rday240 adatrj = adatrj + rdt / rday 241 fjulday = fjulday + rdt / rday 227 242 IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < zprec ) fjulday = REAL(NINT(fjulday),wp) ! avoid truncation error 228 243 IF( ABS(adatrj - REAL(NINT(adatrj ),wp)) < zprec ) adatrj = REAL(NINT(adatrj ),wp) ! avoid truncation error … … 302 317 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 303 318 ! 304 REAL(wp) :: zkt, zndastp 319 REAL(wp) :: zkt, zndastp, zdayfrac, ksecs, ktime 320 INTEGER :: ihour, iminute 305 321 !!---------------------------------------------------------------------- 306 322 … … 327 343 ! define ndastp and adatrj 328 344 IF ( nrstdt == 2 ) THEN 329 ! read the parameters correspond ting to nit000 - 1 (last time step of previous run)345 ! read the parameters corresponding to nit000 - 1 (last time step of previous run) 330 346 CALL iom_get( numror, 'ndastp', zndastp ) 331 347 ndastp = NINT( zndastp ) 332 348 CALL iom_get( numror, 'adatrj', adatrj ) 349 CALL iom_get( numror, 'ntime', ktime ) 350 nn_time0=INT(ktime) 351 ! calculate start time in hours and minutes 352 zdayfrac=adatrj-INT(adatrj) 353 ksecs = NINT(zdayfrac*86400) ! Nearest second to catch rounding errors in adatrj 354 ihour = INT(ksecs/3600) 355 iminute = ksecs/60-ihour*60 356 357 ! Add to nn_time0 358 nhour = nn_time0 / 100 359 nminute = ( nn_time0 - nhour * 100 ) 360 nminute=nminute+iminute 361 362 IF( nminute >= 60 ) THEN 363 nminute=nminute-60 364 nhour=nhour+1 365 ENDIF 366 nhour=nhour+ihour 367 IF( nhour >= 24 ) THEN 368 nhour=nhour-24 369 adatrj=adatrj+1 370 ENDIF 371 nn_time0 = nhour * 100 + nminute 372 adatrj = INT(adatrj) ! adatrj set to integer as nn_time0 updated 333 373 ELSE 334 ! parameters correspondting to nit000 - 1 (as we start the step loop with a call to day) 335 ndastp = ndate0 - 1 ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00 336 adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday 374 ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day) 375 ndastp = ndate0 ! ndate0 read in the namelist in dom_nam 376 nhour = nn_time0 / 100 377 nminute = ( nn_time0 - nhour * 100 ) 378 IF( nhour*3600+nminute*60-ndt05 .lt. 0 ) ndastp=ndastp-1 ! Start hour is specified in the namelist (default 0) 379 adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday 337 380 ! note this is wrong if time step has changed during run 338 381 ENDIF 339 382 ELSE 340 ! parameters correspondting to nit000 - 1 (as we start the step loop with a call to day) 341 ndastp = ndate0 - 1 ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00 342 adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday 383 ! parameters corresponding to nit000 - 1 (as we start the step loop with a call to day) 384 ndastp = ndate0 ! ndate0 read in the namelist in dom_nam 385 nhour = nn_time0 / 100 386 nminute = ( nn_time0 - nhour * 100 ) 387 IF( nhour*3600+nminute*60-ndt05 .lt. 0 ) ndastp=ndastp-1 ! Start hour is specified in the namelist (default 0) 388 adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday 343 389 ENDIF 344 390 IF( ABS(adatrj - REAL(NINT(adatrj),wp)) < 0.1 / rday ) adatrj = REAL(NINT(adatrj),wp) ! avoid truncation error … … 348 394 WRITE(numout,*) ' date ndastp : ', ndastp 349 395 WRITE(numout,*) ' number of elapsed days since the begining of run : ', adatrj 396 WRITE(numout,*) ' nn_time0 : ',nn_time0 350 397 WRITE(numout,*) 351 398 ENDIF … … 363 410 CALL iom_rstput( kt, nitrst, numrow, 'adatrj' , adatrj ) ! number of elapsed days since 364 411 ! ! the begining of the run [s] 412 CALL iom_rstput( kt, nitrst, numrow, 'ntime' , REAL( nn_time0, wp) ) ! time 365 413 ENDIF 366 414 ! -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r6138 r6141 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 … … 32 33 REAL(wp), PUBLIC :: rn_bathy !: depth of flat bottom (active if nn_bathy=0; if =0 depth=jpkm1) 33 34 REAL(wp), PUBLIC :: rn_hmin !: minimum ocean depth (>0) or minimum number of ocean levels (<0) 35 REAL(wp), PUBLIC :: rn_isfhmin !: threshold to discriminate grounded ice to floating ice 34 36 REAL(wp), PUBLIC :: rn_e3zps_min !: miminum thickness for partial steps (meters) 35 37 REAL(wp), PUBLIC :: rn_e3zps_rat !: minimum thickness ration for partial steps 36 38 INTEGER , PUBLIC :: nn_msh !: = 1 create a mesh-mask file 37 INTEGER , PUBLIC :: nn_acc !: = 0/1 use of the acceleration of convergence technique38 39 REAL(wp), PUBLIC :: rn_atfp !: asselin time filter parameter 39 REAL(wp), PUBLIC :: rn_rdt !: time step for the dynamics (and tracer if nacc=0) 40 REAL(wp), PUBLIC :: rn_rdtmin !: minimum time step on tracers 41 REAL(wp), PUBLIC :: rn_rdtmax !: maximum time step on tracers 42 REAL(wp), PUBLIC :: rn_rdth !: depth variation of tracer step 40 REAL(wp), PUBLIC :: rn_rdt !: time step for the dynamics and tracer 43 41 INTEGER , PUBLIC :: nn_closea !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 44 42 INTEGER , PUBLIC :: nn_euler !: =0 start with forward time step or not (=1) 43 LOGICAL , PUBLIC :: ln_iscpl !: coupling with ice sheet 45 44 LOGICAL , PUBLIC :: ln_crs !: Apply grid coarsening to dynamical model output or online passive tracers 46 45 47 46 !! Free surface parameters 48 47 !! ======================= 49 LOGICAL , PUBLIC :: ln_dynspg_exp!: Explicit free surface flag50 LOGICAL , PUBLIC :: ln_dynspg_ts!: Split-Explicit free surface flag48 LOGICAL , PUBLIC :: ln_dynspg_exp !: Explicit free surface flag 49 LOGICAL , PUBLIC :: ln_dynspg_ts !: Split-Explicit free surface flag 51 50 52 51 !! Time splitting parameters … … 61 60 !! Horizontal grid parameters for domhgr 62 61 !! ===================================== 63 INTEGER :: jphgr_msh !: type of horizontal mesh62 INTEGER :: jphgr_msh !: type of horizontal mesh 64 63 ! ! = 0 curvilinear coordinate on the sphere read in coordinate.nc 65 64 ! ! = 1 geographical mesh on the sphere with regular grid-spacing … … 68 67 ! ! = 4 Mercator grid with T/U point at the equator 69 68 70 REAL(wp) :: ppglam0 71 REAL(wp) :: ppgphi0 69 REAL(wp) :: ppglam0 !: longitude of first raw and column T-point (jphgr_msh = 1) 70 REAL(wp) :: ppgphi0 !: latitude of first raw and column T-point (jphgr_msh = 1) 72 71 ! ! 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 72 REAL(wp) :: ppe1_deg !: zonal grid-spacing (degrees) 73 REAL(wp) :: ppe2_deg !: meridional grid-spacing (degrees) 74 REAL(wp) :: ppe1_m !: zonal grid-spacing (degrees) 75 REAL(wp) :: ppe2_m !: meridional grid-spacing (degrees) 77 76 78 77 !! Vertical grid parameter for domzgr 79 78 !! ================================== 80 REAL(wp) :: ppsur 81 REAL(wp) :: ppa0 82 REAL(wp) :: ppa1 83 REAL(wp) :: ppkth 84 REAL(wp) :: ppacr 79 REAL(wp) :: ppsur !: ORCA r4, r2 and r05 coefficients 80 REAL(wp) :: ppa0 !: (default coefficients) 81 REAL(wp) :: ppa1 !: 82 REAL(wp) :: ppkth !: 83 REAL(wp) :: ppacr !: 85 84 ! 86 85 ! If both ppa0 ppa1 and ppsur are specified to 0, then 87 86 ! they are computed from ppdzmin, pphmax , ppkth, ppacr in dom_zgr 88 REAL(wp) :: ppdzmin 89 REAL(wp) :: pphmax 87 REAL(wp) :: ppdzmin !: Minimum vertical spacing 88 REAL(wp) :: pphmax !: Maximum depth 90 89 ! 91 LOGICAL :: ldbletanh 92 REAL(wp) :: ppa2 93 REAL(wp) :: ppkth2 94 REAL(wp) :: ppacr2 90 LOGICAL :: ldbletanh !: Use/do not use double tanf function for vertical coordinates 91 REAL(wp) :: ppa2 !: Double tanh function parameters 92 REAL(wp) :: ppkth2 !: 93 REAL(wp) :: ppacr2 !: 95 94 96 95 ! !! old non-DOCTOR names still used in the model … … 99 98 REAL(wp), PUBLIC :: e3zps_rat !: minimum thickness ration for partial steps 100 99 INTEGER , PUBLIC :: nmsh !: = 1 create a mesh-mask file 101 INTEGER , PUBLIC :: nacc !: = 0/1 use of the acceleration of convergence technique102 100 REAL(wp), PUBLIC :: atfp !: asselin time filter parameter 103 REAL(wp), PUBLIC :: rdt !: time step for the dynamics (and tracer if nacc=0) 104 REAL(wp), PUBLIC :: rdtmin !: minimum time step on tracers 105 REAL(wp), PUBLIC :: rdtmax !: maximum time step on tracers 106 REAL(wp), PUBLIC :: rdth !: depth variation of tracer step 101 REAL(wp), PUBLIC :: rdt !: time step for the dynamics and tracer 107 102 108 103 ! !!! associated variables 109 104 INTEGER , PUBLIC :: neuler !: restart euler forward option (0=Euler) 110 105 REAL(wp), PUBLIC :: atfp1 !: asselin time filter coeff. (atfp1= 1-2*atfp) 111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: rdttra !: vertical profile of tracer time step 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: r2dtra !: = 2*rdttra except at nit000 (=rdttra) if neuler=0 106 REAL(wp), PUBLIC :: r2dt !: = 2*rdt except at nit000 (=rdt) if neuler=0 113 107 114 108 !!---------------------------------------------------------------------- … … 177 171 !! vertical coordinate and scale factors 178 172 !! --------------------------------------------------------------------- 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) 173 ! !!* Namelist namzgr : vertical coordinate * 174 LOGICAL, PUBLIC :: ln_zco !: z-coordinate - full step 175 LOGICAL, PUBLIC :: ln_zps !: z-coordinate - partial step 176 LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate 177 LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF 178 LOGICAL, PUBLIC :: ln_linssh !: variable grid flag 179 180 ! ! ref. ! before ! now ! after ! 181 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 , e3t_b , e3t_n , e3t_a !: t- vert. scale factor [m] 182 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_0 , e3u_b , e3u_n , e3u_a !: u- vert. scale factor [m] 183 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 , e3v_b , e3v_n , e3v_a !: v- vert. scale factor [m] 184 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f_0 , e3f_n !: f- vert. scale factor [m] 185 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 , e3w_b , e3w_n !: w- vert. scale factor [m] 186 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_0 , e3uw_b , e3uw_n !: uw-vert. scale factor [m] 187 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 , e3vw_b , e3vw_n !: vw-vert. scale factor [m] 188 189 ! ! ref. ! before ! now ! 190 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 , gdept_b , gdept_n !: t- depth [m] 191 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdepw_0 , gdepw_b , gdepw_n !: w- depth [m] 192 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0 , gde3w_n !: w- depth (sum of e3w) [m] 193 194 ! ! ref. ! before ! now ! after ! 195 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 , ht_n !: t-depth [m] 196 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hu_b , hu_n , hu_a !: u-depth [m] 197 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0 , hv_b , hv_n , hv_a !: u-depth [m] 198 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_b , r1_hu_n , r1_hu_a !: inverse of u-depth [1/m] 199 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_b , r1_hv_n , r1_hv_a !: inverse of v-depth [1/m] 200 223 201 224 202 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) 225 203 INTEGER, PUBLIC :: nlb10 !: shallowest W level Bellow ~10m (nla10 + 1) 226 204 227 !! z-coordinate with full steps (also used in the other cases as reference z-coordinate)205 !! 1D reference vertical coordinate 228 206 !! =-----------------====------ 229 207 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m) … … 231 209 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3tp , e3wp !: ocean bottom level thickness at T and W points 232 210 211 !!gm This should be removed from here.... ==>>> only used in domzgr at initialization phase 233 212 !! s-coordinate and hybrid z-s-coordinate 234 213 !! =----------------======--------------- … … 244 223 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing t--u points (m) 245 224 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio 225 !!gm end 246 226 247 227 !!---------------------------------------------------------------------- … … 251 231 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level 252 232 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 254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i , umask_i, vmask_i, fmask_i!: interior domain T-point mask255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function233 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 234 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 235 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_h !: internal domain T-point mask (Figure 8.5 NEMO book) 256 236 257 237 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfdep !: top first ocean level (ISF) 258 238 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: first wet T-, U-, V-, F- ocean level (ISF) 259 239 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep !: Iceshelf draft (ISF) 260 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask !: surface domain T-point mask 261 240 241 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts 262 242 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 263 243 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask !: land/ocean mask at WT-, WU- and WV-pts … … 271 251 INTEGER , PUBLIC :: nmonth !: current month 272 252 INTEGER , PUBLIC :: nday !: current day of the month 253 INTEGER , PUBLIC :: nhour !: current hour 254 INTEGER , PUBLIC :: nminute !: current minute 273 255 INTEGER , PUBLIC :: ndastp !: time step date in yyyymmdd format 274 256 INTEGER , PUBLIC :: nday_year !: current day counted from jan 1st of the current year … … 307 289 !!---------------------------------------------------------------------- 308 290 !! NEMO/OPA 4.0 , NEMO Consortium (2011) 309 !! $Id$ 291 !! $Id$ 310 292 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 311 293 !!---------------------------------------------------------------------- … … 331 313 ierr(:) = 0 332 314 ! 333 ALLOCATE( rdttra(jpk), r2dtra(jpk),mig(jpi), mjg(jpj), nfiimpp(jpni,jpnj), &315 ALLOCATE( mig(jpi), mjg(jpj), nfiimpp(jpni,jpnj), & 334 316 & nfipproc(jpni,jpnj), nfilcit(jpni,jpnj), STAT=ierr(1) ) 335 317 ! … … 352 334 & ff (jpi,jpj) , STAT=ierr(3) ) 353 335 ! 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) ) 336 ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) , & 337 & gdept_b(jpi,jpj,jpk) , gdepw_b(jpi,jpj,jpk) , & 338 & gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) ) 339 ! 340 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) , & 341 & e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) , e3w_b(jpi,jpj,jpk) , & 342 & 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) , & 343 & e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) , & 344 ! ! 345 & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , & 346 & e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) , & 347 & e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , STAT=ierr(5) ) 348 ! 349 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , & 350 & hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) , & 351 & ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) , & 352 & hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6) ) 353 ! 374 354 ! 375 355 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , & … … 384 364 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) ) 385 365 386 ALLOCATE( mbathy(jpi,jpj) , bathy (jpi,jpj) ,&387 & tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), &388 & bmask (jpi,jpj) ,&389 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) )366 ALLOCATE( mbathy(jpi,jpj) , bathy (jpi,jpj) , & 367 & tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 368 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 369 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) ) 390 370 391 371 ! (ISF) Allocation of basic array 392 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), 393 & 394 & mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) )372 ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj), & 373 & mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) , & 374 & mikf(jpi,jpj), STAT=ierr(10) ) 395 375 396 376 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), & … … 406 386 !!====================================================================== 407 387 END MODULE dom_oce 408 -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5870 r6141 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(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 122 z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 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 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file 143 IF( nmsh /= 0 .AND. ln_iscpl .AND. .NOT. ln_rstart ) CALL dom_wri ! Create a domain file 144 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 120 145 ! 121 146 IF( nn_timing == 1 ) CALL timing_stop('dom_init') … … 135 160 !!---------------------------------------------------------------------- 136 161 USE ioipsl 137 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, & 138 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl, & 139 & nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , & 140 & nn_write, ln_dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler 141 NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, & 142 & nn_acc , rn_atfp , rn_rdt , rn_rdtmin , & 143 & rn_rdtmax, rn_rdth , nn_closea , ln_crs, & 144 & jphgr_msh, & 145 & ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, & 146 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 162 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, & 163 nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl , & 164 & nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate , & 165 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, nn_euler , & 166 & ln_cfmeta, ln_iscpl 167 NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, & 168 & rn_atfp , rn_rdt , nn_closea , ln_crs , jphgr_msh , & 169 & ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, & 170 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 147 171 & ppa2, ppkth2, ppacr2 148 172 #if defined key_netcdf4 … … 178 202 WRITE(numout,*) ' number of the last time step nn_itend = ', nn_itend 179 203 WRITE(numout,*) ' initial calendar date aammjj nn_date0 = ', nn_date0 204 WRITE(numout,*) ' initial time of day in hhmm nn_time0 = ', nn_time0 180 205 WRITE(numout,*) ' leap year calendar (0/1) nn_leapy = ', nn_leapy 181 206 WRITE(numout,*) ' initial state output nn_istate = ', nn_istate … … 186 211 ENDIF 187 212 WRITE(numout,*) ' frequency of output file nn_write = ', nn_write 188 WRITE(numout,*) ' multi file dimgout ln_dimgnnn = ', ln_dimgnnn189 213 WRITE(numout,*) ' mask land points ln_mskland = ', ln_mskland 190 214 WRITE(numout,*) ' additional CF standard metadata ln_cfmeta = ', ln_cfmeta 191 215 WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber 192 216 WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz 217 WRITE(numout,*) ' IS coupling at the restart step ln_iscpl = ', ln_iscpl 193 218 ENDIF 194 219 … … 258 283 WRITE(numout,*) ' min depth of the ocean (>0) or rn_hmin = ', rn_hmin 259 284 WRITE(numout,*) ' min number of ocean level (<0) ' 285 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' (m)' 260 286 WRITE(numout,*) ' minimum thickness of partial rn_e3zps_min = ', rn_e3zps_min, ' (m)' 261 287 WRITE(numout,*) ' step level rn_e3zps_rat = ', rn_e3zps_rat … … 267 293 WRITE(numout,*) ' ocean time step rn_rdt = ', rn_rdt 268 294 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 269 WRITE(numout,*) ' acceleration of converge nn_acc = ', nn_acc270 WRITE(numout,*) ' nn_acc=1: surface tracer rdt rn_rdtmin = ', rn_rdtmin271 WRITE(numout,*) ' bottom tracer rdt rdtmax = ', rn_rdtmax272 WRITE(numout,*) ' depth of transition rn_rdth = ', rn_rdth273 295 WRITE(numout,*) ' suppression of closed seas (=0) nn_closea = ', nn_closea 274 296 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs … … 297 319 e3zps_rat = rn_e3zps_rat 298 320 nmsh = nn_msh 299 nacc = nn_acc300 321 atfp = rn_atfp 301 322 rdt = rn_rdt 302 rdtmin = rn_rdtmin303 rdtmax = rn_rdtmin304 rdth = rn_rdth305 323 306 324 #if defined key_netcdf4 … … 403 421 INTEGER :: ji, jj, jk 404 422 REAL(wp) :: zrxmax 405 REAL(wp), DIMENSION(4) :: zr1423 REAL(wp), DIMENSION(4) :: zr1 406 424 !!---------------------------------------------------------------------- 407 425 rx1(:,:) = 0._wp … … 412 430 DO jj = 2, jpjm1 413 431 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)432 zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) & 433 & +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) & 434 & / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) & 435 & -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk) 436 zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) & 437 & +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) & 438 & / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) & 439 & -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk) 440 zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) & 441 & +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) & 442 & / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) & 443 & -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk) 444 zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) & 445 & +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) & 446 & / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) & 447 & -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk) 448 zrxmax = MAXVAL( zr1(1:4) ) 449 rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax ) 432 450 END DO 433 451 END DO -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90
r4667 r6141 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_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r5870 r6141 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 … … 506 506 CALL iom_close( inum ) 507 507 508 !!gm THIS is TO BE REMOVED !!!!!!!509 510 ! need to be define for the extended grid south of -80S511 ! some point are undefined but you need to have e1 and e2 .NE. 0512 WHERE (e1t==0.0_wp)513 e1t=1.0e2514 END WHERE515 WHERE (e1v==0.0_wp)516 e1v=1.0e2517 END WHERE518 WHERE (e1u==0.0_wp)519 e1u=1.0e2520 END WHERE521 WHERE (e1f==0.0_wp)522 e1f=1.0e2523 END WHERE524 WHERE (e2t==0.0_wp)525 e2t=1.0e2526 END WHERE527 WHERE (e2v==0.0_wp)528 e2v=1.0e2529 END WHERE530 WHERE (e2u==0.0_wp)531 e2u=1.0e2532 END WHERE533 WHERE (e2f==0.0_wp)534 e2f=1.0e2535 END WHERE536 !!gm end537 538 508 END SUBROUTINE hgr_read 539 509 -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r6138 r6141 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 !!---------------------------------------------------------------------- … … 183 174 ! -------------------- 184 175 tmask_i(:,:) = ssmask(:,:) ! (ISH) tmask_i = 1 even on the ice shelf 176 177 tmask_h(:,:) = 1._wp ! 0 on the halo and 1 elsewhere 185 178 iif = jpreci ! ??? 186 179 iil = nlci - jpreci + 1 … … 188 181 ijl = nlcj - jprecj + 1 189 182 190 tmask_ i( 1 :iif, : ) = 0._wp ! first columns191 tmask_ i(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns)192 tmask_ i( : , 1 :ijf) = 0._wp ! first rows193 tmask_ i( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows)183 tmask_h( 1 :iif, : ) = 0._wp ! first columns 184 tmask_h(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns) 185 tmask_h( : , 1 :ijf) = 0._wp ! first rows 186 tmask_h( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows) 194 187 195 188 ! north fold mask … … 202 195 IF( mjg(nlej) == jpjglo ) THEN ! only half of the nlcj-1 row 203 196 DO ji = iif+1, iil-1 204 tmask_ i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))197 tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji)) 205 198 END DO 206 199 ENDIF 207 200 ENDIF 201 202 tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:) 203 208 204 IF( jperio == 5 .OR. jperio == 6 ) THEN ! F-point pivot 209 205 tpol( 1 :jpiglo) = 0._wp … … 225 221 END DO 226 222 END DO 227 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point223 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point 228 224 DO jj = 1, jpjm1 229 225 DO ji = 1, fs_jpim1 ! vector loop 230 umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:)))231 vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:)))226 ssumask(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:))) 227 ssvmask(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:))) 232 228 END DO 233 229 DO ji = 1, jpim1 ! NO vector opt. 234 fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) &230 ssfmask(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 235 231 & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 236 232 END DO 237 233 END DO 238 CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions239 CALL lbc_lnk( vmask , 'V', 1._wp )240 CALL lbc_lnk( fmask , 'F', 1._wp )241 CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions242 CALL lbc_lnk( vmask_i, 'V', 1._wp )243 CALL lbc_lnk( fmask_i, 'F', 1._wp )234 CALL lbc_lnk( umask , 'U', 1._wp ) ! Lateral boundary conditions 235 CALL lbc_lnk( vmask , 'V', 1._wp ) 236 CALL lbc_lnk( fmask , 'F', 1._wp ) 237 CALL lbc_lnk( ssumask, 'U', 1._wp ) ! Lateral boundary conditions 238 CALL lbc_lnk( ssvmask, 'V', 1._wp ) 239 CALL lbc_lnk( ssfmask, 'F', 1._wp ) 244 240 245 241 ! 3. Ocean/land mask at wu-, wv- and w points … … 253 249 wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 254 250 END DO 255 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 251 302 252 ! Lateral boundary conditions on velocity (modify fmask) … … 399 349 ! 400 350 CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask 401 351 ! 402 352 ! 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 353 ! 449 354 CALL wrk_dealloc( jpi, jpj, imsk ) -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
r5870 r6141 28 28 CONTAINS 29 29 30 SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid )30 SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk ) 31 31 !!---------------------------------------------------------------------- 32 32 !! *** ROUTINE dom_ngb *** … … 39 39 REAL(wp) , INTENT(in ) :: plon, plat ! longitude,latitude of the point 40 40 INTEGER , INTENT( out) :: kii, kjj ! i-,j-index of the closes grid point 41 INTEGER , INTENT(in ), OPTIONAL :: kkk ! k-index of the mask level used 41 42 CHARACTER(len=1), INTENT(in ) :: cdgrid ! grid name 'T', 'U', 'V', 'W' 42 43 ! 44 INTEGER :: ik ! working level 43 45 INTEGER , DIMENSION(2) :: iloc 44 46 REAL(wp) :: zlon, zmini … … 51 53 ! 52 54 zmask(:,:) = 0._wp 55 ik = 1 56 IF ( PRESENT(kkk) ) ik=kkk 53 57 SELECT CASE( cdgrid ) 54 CASE( 'U' ) ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej, 1)55 CASE( 'V' ) ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej, 1)56 CASE( 'F' ) ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej, 1)57 CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej, 1)58 CASE( 'U' ) ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik) 59 CASE( 'V' ) ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik) 60 CASE( 'F' ) ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik) 61 CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik) 58 62 END SELECT 59 63 60 zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360 61 zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360 62 IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270 63 IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180 64 IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN 65 zlon = MOD( plon + 720., 360. ) ! plon between 0 and 360 66 zglam(:,:) = MOD( zglam(:,:) + 720., 360. ) ! glam between 0 and 360 67 IF( zlon > 270. ) zlon = zlon - 360. ! zlon between -90 and 270 68 IF( zlon < 90. ) WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360. ! glam between -180 and 180 69 zglam(:,:) = zglam(:,:) - zlon 70 ELSE 71 zglam(:,:) = zglam(:,:) - plon 72 END IF 64 73 65 zglam(:,:) = zglam(:,:) - zlon66 74 zgphi(:,:) = zgphi(:,:) - plat 67 75 zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:) -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domstp.F90
r4292 r6141 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) … … 41 39 !! filter parameter read in namelist 42 40 !! - Model time step: 43 !! nacc = 0 : synchronous time intergration. 44 !! There is one time step only, defined by: rdt, rdttra(k)=rdt 45 !! nacc = 1 : accelerating the convergence. There is 2 different 46 !! time steps for dynamics and tracers: 47 !! rdt : dynamical part 48 !! rdttra(k): temperature and salinity 49 !! The tracer time step is a function of vertical level. the model 50 !! reference time step ( i.e. for wind stress, surface heat and 51 !! salt fluxes) is the surface tracer time step is rdttra(1). 52 !! N.B. depth dependent acceleration of convergence is not im- 53 !! plemented for s-coordinate. 41 !! synchronous time intergration. 42 !! There is one time step only, defined by: rdt for dynamics and 43 !! tracer,wind stress, surface heat and salt fluxes 54 44 !! 55 !! ** Action : - rdttra : vertical profile of tracer time step45 !! ** Action : [REMOVED - rdttra: vertical profile of tracer time step] 56 46 !! - atfp1 : = 1 - 2*atfp 57 47 !! … … 72 62 atfp1 = 1. - 2. * atfp 73 63 74 SELECT CASE ( nacc ) 64 IF(lwp) WRITE(numout,*)' synchronous time stepping' 65 IF(lwp) WRITE(numout,*)' dynamics and tracer time step = ', rdt/3600., ' hours' 75 66 76 CASE ( 0 ) ! Synchronous time stepping77 IF(lwp) WRITE(numout,*)' synchronous time stepping'78 IF(lwp) WRITE(numout,*)' dynamics and tracer time step = ', rdt/3600., ' hours'79 80 rdttra(:) = rdt81 82 CASE ( 1 ) ! Accelerating the convergence83 IF(lwp) WRITE(numout,*) ' no tracer damping in the turbocline'84 IF(lwp) WRITE(numout,*)' accelerating the convergence'85 IF(lwp) WRITE(numout,*)' dynamics time step = ', rdt/3600., ' hours'86 IF( ln_sco .AND. rdtmin /= rdtmax .AND. lk_vvl ) &87 & CALL ctl_stop ( ' depth dependent acceleration of convergence not implemented in s-coordinates &88 & nor in variable volume' )89 IF(lwp) WRITE(numout,*)' tracers time step : dt (hours) level'90 91 DO jk = 1, jpk92 IF( gdept_1d(jk) <= rdth ) rdttra(jk) = rdtmin93 IF( gdept_1d(jk) > rdth ) THEN94 rdttra(jk) = rdtmin + ( rdtmax - rdtmin ) &95 * ( EXP( ( gdept_1d(jk ) - rdth ) / rdth ) - 1. ) &96 / ( EXP( ( gdept_1d(jpk) - rdth ) / rdth ) - 1. )97 ENDIF98 IF(lwp) WRITE(numout,"(36x,f5.2,5x,i3)") rdttra(jk)/3600., jk99 END DO100 101 CASE DEFAULT ! E R R O R102 103 WRITE(ctmp1,*) ' nacc value e r r o r, nacc= ',nacc104 CALL ctl_stop( ctmp1 )105 106 END SELECT107 67 108 68 END SUBROUTINE dom_stp -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5870 r6141 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 24 25 USE wet_dry ! wetting and drying 26 USE restart ! ocean restart 27 ! 25 28 USE in_out_manager ! I/O manager 26 29 USE iom ! I/O manager library 27 USE restart ! ocean restart28 30 USE lib_mpp ! distributed memory computing library 29 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 61 63 62 64 !! * Substitutions 63 # include "domzgr_substitute.h90"64 65 # include "vectopt_loop_substitute.h90" 65 66 !!---------------------------------------------------------------------- 66 !! NEMO/OPA 3. 3 , NEMO-Consortium (2010)67 !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 67 68 !! $Id$ 68 69 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 69 70 !!---------------------------------------------------------------------- 70 71 71 CONTAINS 72 72 … … 75 75 !! *** FUNCTION dom_vvl_alloc *** 76 76 !!---------------------------------------------------------------------- 77 IF( ln_vvl_zstar ) dom_vvl_alloc = 077 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 78 78 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 79 79 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & … … 82 82 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 83 83 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 84 un_td = 0. 0_wp85 vn_td = 0. 0_wp84 un_td = 0._wp 85 vn_td = 0._wp 86 86 ENDIF 87 87 IF( ln_vvl_ztilde ) THEN … … 90 90 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 91 91 ENDIF 92 92 ! 93 93 END FUNCTION dom_vvl_alloc 94 94 … … 104 104 !! - interpolate scale factors 105 105 !! 106 !! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b)107 !! - Regrid: fse3(u/v)_n108 !! fse3(u/v)_b109 !! fse3w_n110 !! fse3(u/v)w_b111 !! fse3(u/v)w_n112 !! fsdept_n, fsdepw_n and fsde3w_n106 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 107 !! - Regrid: e3(u/v)_n 108 !! e3(u/v)_b 109 !! e3w_n 110 !! e3(u/v)w_b 111 !! e3(u/v)w_n 112 !! gdept_n, gdepw_n and gde3w_n 113 113 !! - h(t/u/v)_0 114 114 !! - frq_rst_e3t and frq_rst_hdv … … 116 116 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 117 117 !!---------------------------------------------------------------------- 118 USE phycst, ONLY : rpi, rsmall, rad 119 !! * Local declarations 120 INTEGER :: ji,jj,jk 118 INTEGER :: ji, jj, jk 121 119 INTEGER :: ii0, ii1, ij0, ij1 122 120 REAL(wp):: zcoef 123 121 !!---------------------------------------------------------------------- 124 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 125 122 ! 123 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 124 ! 126 125 IF(lwp) WRITE(numout,*) 127 126 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 128 127 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 129 130 ! choose vertical coordinate (z_star, z_tilde or layer) 131 ! ========================== 132 CALL dom_vvl_ctl 133 134 ! Allocate module arrays 135 ! ====================== 128 ! 129 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 130 ! 131 ! ! Allocate module arrays 136 132 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 137 138 ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk)) 139 ! ============================================================================ 133 ! 134 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 140 135 CALL dom_vvl_rst( nit000, 'READ' ) 141 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 142 143 ! Reconstruction of all vertical scale factors at now and before time steps 144 ! ============================================================================= 145 ! Horizontal scale factor interpolations 146 ! -------------------------------------- 147 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 148 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 149 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 150 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 151 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 152 ! Vertical scale factor interpolations 153 ! ------------------------------------ 154 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 155 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 156 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 157 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 158 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 159 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 160 ! t- and w- points depth 161 ! ---------------------- 162 ! set the isf depth as it is in the initial step 163 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 164 fsdepw_n(:,:,1) = 0.0_wp 165 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 166 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 167 fsdepw_b(:,:,1) = 0.0_wp 168 169 DO jk = 2, jpk 136 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 137 ! 138 ! !== Set of all other vertical scale factors ==! (now and before) 139 ! ! Horizontal interpolation of e3t 140 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U 141 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 142 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V 143 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 144 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F 145 ! ! Vertical interpolation of e3t,u,v 146 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W 147 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' ) 148 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW 149 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 150 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW 151 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 152 ! 153 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 154 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 155 gdepw_n(:,:,1) = 0.0_wp 156 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 157 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 158 gdepw_b(:,:,1) = 0.0_wp 159 DO jk = 2, jpk ! vertical sum 170 160 DO jj = 1,jpj 171 161 DO ji = 1,jpi 172 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 173 ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 174 ! 0.5 where jk = mikt 175 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 176 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 177 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 178 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 179 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 180 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 181 fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) & 182 & + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)) 162 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 163 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 164 ! ! 0.5 where jk = mikt 165 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 166 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 167 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 168 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 169 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 170 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 171 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 172 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 173 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 183 174 END DO 184 175 END DO 185 176 END DO 186 187 ! Before depth and Inverse of the local depth of the water column at u- and v- points 188 ! ----------------------------------------------------------------------------------- 189 hu_b(:,:) = 0. 190 hv_b(:,:) = 0. 191 DO jk = 1, jpkm1 192 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 193 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 177 ! 178 ! !== thickness of the water column !! (ocean portion only) 179 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 180 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 181 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 182 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 183 hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1) 184 DO jk = 2, jpkm1 185 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 186 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 187 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 188 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 189 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 194 190 END DO 195 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 196 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 197 198 ! Restoring frequencies for z_tilde coordinate 199 ! ============================================ 191 ! 192 ! !== inverse of water column thickness ==! (u- and v- points) 193 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 194 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 195 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 196 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 197 198 ! !== z_tilde coordinate case ==! (Restoring frequencies) 200 199 IF( ln_vvl_ztilde ) THEN 201 ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 202 frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 203 frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 204 IF( ln_vvl_ztilde_as_zstar ) THEN 205 ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 206 frq_rst_e3t(:,:) = 0.0_wp 207 frq_rst_hdv(:,:) = 1.0_wp / rdt 200 !!gm : idea: add here a READ in a file of custumized restoring frequency 201 ! ! Values in days provided via the namelist 202 ! ! use rsmall to avoid possible division by zero errors with faulty settings 203 frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 204 frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 205 ! 206 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 207 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 208 frq_rst_hdv(:,:) = 1._wp / rdt 208 209 ENDIF 209 IF ( ln_vvl_zstar_at_eqtor ) THEN 210 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 210 211 DO jj = 1, jpj 211 212 DO ji = 1, jpi 213 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 212 214 IF( ABS(gphit(ji,jj)) >= 6.) THEN 213 215 ! values outside the equatorial band and transition zone (ztilde) 214 216 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 215 217 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 216 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 218 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 217 219 ! values inside the equatorial band (ztilde as zstar) 218 220 frq_rst_e3t(ji,jj) = 0.0_wp 219 221 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 220 ELSE 221 ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)222 ELSE ! transition band (2.5 to 6 degrees N/S) 223 ! ! (linearly transition from z-tilde to z-star) 222 224 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 223 225 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & … … 230 232 END DO 231 233 END DO 232 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN 233 ii0 = 103 ; ii1 = 111 ! Suppress ztilde in the Foxe Basin for ORCA2234 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 235 ii0 = 103 ; ii1 = 111 234 236 ij0 = 128 ; ij1 = 135 ; 235 237 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp … … 238 240 ENDIF 239 241 ENDIF 240 242 ! 241 243 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 242 244 ! 243 245 END SUBROUTINE dom_vvl_init 244 246 … … 262 264 !! - tilde_e3t_a: after increment of vertical scale factor 263 265 !! in z_tilde case 264 !! - fse3(t/u/v)_a266 !! - e3(t/u/v)_a 265 267 !! 266 268 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 267 269 !!---------------------------------------------------------------------- 268 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 269 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 270 !! * Arguments 271 INTEGER, INTENT( in ) :: kt ! time step 272 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 273 !! * Local declarations 274 INTEGER :: ji, jj, jk ! dummy loop indices 275 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 276 REAL(wp) :: z2dt ! temporary scalars 277 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 278 LOGICAL :: ll_do_bclinic ! temporary logical 279 !!---------------------------------------------------------------------- 280 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 281 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 282 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 283 284 IF(kt == nit000) THEN 270 INTEGER, INTENT( in ) :: kt ! time step 271 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 272 ! 273 INTEGER :: ji, jj, jk ! dummy loop indices 274 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 275 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars 276 LOGICAL :: ll_do_bclinic ! local logical 277 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 278 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 279 !!---------------------------------------------------------------------- 280 ! 281 IF( ln_linssh ) RETURN ! No calculation in linear free surface 282 ! 283 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 284 ! 285 CALL wrk_alloc( jpi,jpj,zht, z_scale, zwu, zwv, zhdiv ) 286 CALL wrk_alloc( jpi,jpj,jpk, ze3t ) 287 288 IF( kt == nit000 ) THEN 285 289 IF(lwp) WRITE(numout,*) 286 290 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' … … 290 294 ll_do_bclinic = .TRUE. 291 295 IF( PRESENT(kcall) ) THEN 292 IF ( kcall == 2 .AND. ln_vvl_ztilde )ll_do_bclinic = .FALSE.296 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 293 297 ENDIF 294 298 … … 296 300 ! After acale factors at t-points ! 297 301 ! ******************************* ! 298 299 302 ! ! --------------------------------------------- ! 300 303 ! ! z_star coordinate and barotropic z-tilde part ! 301 304 ! ! --------------------------------------------- ! 302 305 ! 303 306 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 304 307 DO jk = 1, jpkm1 305 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)306 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)308 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 309 e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 307 310 END DO 308 311 ! 309 312 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 310 313 ! ! ------baroclinic part------ ! 311 312 314 ! I - initialization 313 315 ! ================== … … 315 317 ! 1 - barotropic divergence 316 318 ! ------------------------- 317 zhdiv(:,:) = 0. 318 zht(:,:) = 0. 319 zhdiv(:,:) = 0._wp 320 zht(:,:) = 0._wp 319 321 DO jk = 1, jpkm1 320 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)321 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)322 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 323 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 322 324 END DO 323 325 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 326 328 ! -------------------------------------------------- 327 329 IF( ln_vvl_ztilde ) THEN 328 IF( kt .GT.nit000 ) THEN330 IF( kt > nit000 ) THEN 329 331 DO jk = 1, jpkm1 330 332 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 331 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )333 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 332 334 END DO 333 335 ENDIF 334 END 336 ENDIF 335 337 336 338 ! II - after z_tilde increments of vertical scale factors 337 339 ! ======================================================= 338 tilde_e3t_a(:,:,:) = 0. 0_wp ! tilde_e3t_a used to store tendency terms340 tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms 339 341 340 342 ! 1 - High frequency divergence term … … 342 344 IF( ln_vvl_ztilde ) THEN ! z_tilde case 343 345 DO jk = 1, jpkm1 344 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )346 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 345 347 END DO 346 348 ELSE ! layer case 347 349 DO jk = 1, jpkm1 348 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)349 END DO 350 END 350 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 351 END DO 352 ENDIF 351 353 352 354 ! 2 - Restoring term (z-tilde case only) … … 356 358 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 357 359 END DO 358 END 360 ENDIF 359 361 360 362 ! 3 - Thickness diffusion term 361 363 ! ---------------------------- 362 zwu(:,:) = 0.0_wp 363 zwv(:,:) = 0.0_wp 364 ! a - first derivative: diffusive fluxes 365 DO jk = 1, jpkm1 364 zwu(:,:) = 0._wp 365 zwv(:,:) = 0._wp 366 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 366 367 DO jj = 1, jpjm1 367 368 DO ji = 1, fs_jpim1 ! vector opt. … … 375 376 END DO 376 377 END DO 377 ! b - correction for last oceanic u-v points 378 DO jj = 1, jpj 378 DO jj = 1, jpj ! b - correction for last oceanic u-v points 379 379 DO ji = 1, jpi 380 380 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) … … 382 382 END DO 383 383 END DO 384 ! c - second derivative: divergence of diffusive fluxes 385 DO jk = 1, jpkm1 384 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 386 385 DO jj = 2, jpjm1 387 386 DO ji = fs_2, fs_jpim1 ! vector opt. … … 392 391 END DO 393 392 END DO 394 ! d - thickness diffusion transport: boundary conditions395 ! (stored for tracer advction and continuity equation)393 ! ! d - thickness diffusion transport: boundary conditions 394 ! (stored for tracer advction and continuity equation) 396 395 CALL lbc_lnk( un_td , 'U' , -1._wp) 397 396 CALL lbc_lnk( vn_td , 'V' , -1._wp) … … 411 410 ! Maximum deformation control 412 411 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 413 ze3t(:,:,jpk) = 0. 0_wp412 ze3t(:,:,jpk) = 0._wp 414 413 DO jk = 1, jpkm1 415 414 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) … … 420 419 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 421 420 ! - ML - test: for the moment, stop simulation for too large e3_t variations 422 IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT.- rn_zdef_max ) ) THEN421 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 423 422 IF( lk_mpp ) THEN 424 423 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) … … 463 462 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 464 463 DO jk = 1, jpkm1 465 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)464 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 466 465 END DO 467 466 … … 471 470 ! ! ---baroclinic part--------- ! 472 471 DO jk = 1, jpkm1 473 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)472 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 474 473 END DO 475 474 ENDIF … … 486 485 zht(:,:) = 0.0_wp 487 486 DO jk = 1, jpkm1 488 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)487 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 489 488 END DO 490 489 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 491 490 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 492 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM( fse3t_n))) =', z_tmax491 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 493 492 ! 494 493 zht(:,:) = 0.0_wp 495 494 DO jk = 1, jpkm1 496 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)495 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 497 496 END DO 498 497 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 499 498 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 500 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM( fse3t_a))) =', z_tmax499 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 501 500 ! 502 501 zht(:,:) = 0.0_wp 503 502 DO jk = 1, jpkm1 504 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)503 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 505 504 END DO 506 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 507 506 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 508 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM( fse3t_b))) =', z_tmax507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 509 508 ! 510 509 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) … … 525 524 ! *********************************** ! 526 525 527 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )528 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )526 CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 527 CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 529 528 530 529 ! *********************************** ! … … 532 531 ! *********************************** ! 533 532 534 hu_a(:,:) = 0._wp ! Ocean depth at U-points535 hv_a(:,:) = 0._wp ! Ocean depth at V-points536 DO jk = 1, jpkm1537 hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)538 hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)533 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 534 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 535 DO jk = 2, jpkm1 536 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 537 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 539 538 END DO 540 539 ! ! Inverse of the local depth 541 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 542 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 543 544 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 545 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 546 540 !!gm BUG ? don't understand the use of umask_i here ..... 541 r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 542 r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 543 ! 544 CALL wrk_dealloc( jpi,jpj, zht, z_scale, zwu, zwv, zhdiv ) 545 CALL wrk_dealloc( jpi,jpj,jpk, ze3t ) 546 ! 547 547 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt') 548 548 ! 549 549 END SUBROUTINE dom_vvl_sf_nxt 550 550 … … 562 562 !! - recompute depths and water height fields 563 563 !! 564 !! ** Action : - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step564 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step 565 565 !! - Recompute: 566 !! fse3(u/v)_b567 !! fse3w_n568 !! fse3(u/v)w_b569 !! fse3(u/v)w_n570 !! fsdept_n, fsdepw_n and fsde3w_n566 !! e3(u/v)_b 567 !! e3w_n 568 !! e3(u/v)w_b 569 !! e3(u/v)w_n 570 !! gdept_n, gdepw_n and gde3w_n 571 571 !! h(u/v) and h(u/v)r 572 572 !! … … 574 574 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 575 575 !!---------------------------------------------------------------------- 576 !! * Arguments 577 INTEGER, INTENT( in ) :: kt ! time step 578 !! * Local declarations 579 INTEGER :: ji,jj,jk ! dummy loop indices 580 REAL(wp) :: zcoef 581 !!---------------------------------------------------------------------- 582 576 INTEGER, INTENT( in ) :: kt ! time step 577 ! 578 INTEGER :: ji, jj, jk ! dummy loop indices 579 REAL(wp) :: zcoef ! local scalar 580 !!---------------------------------------------------------------------- 581 ! 582 IF( ln_linssh ) RETURN ! No calculation in linear free surface 583 ! 583 584 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 584 585 ! … … 591 592 ! Time filter and swap of scale factors 592 593 ! ===================================== 593 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.594 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 594 595 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 595 596 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 601 602 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 602 603 ENDIF 603 fsdept_b(:,:,:) = fsdept_n(:,:,:)604 fsdepw_b(:,:,:) = fsdepw_n(:,:,:)605 606 fse3t_n(:,:,:) = fse3t_a(:,:,:)607 fse3u_n(:,:,:) = fse3u_a(:,:,:)608 fse3v_n(:,:,:) = fse3v_a(:,:,:)604 gdept_b(:,:,:) = gdept_n(:,:,:) 605 gdepw_b(:,:,:) = gdepw_n(:,:,:) 606 607 e3t_n(:,:,:) = e3t_a(:,:,:) 608 e3u_n(:,:,:) = e3u_a(:,:,:) 609 e3v_n(:,:,:) = e3v_a(:,:,:) 609 610 610 611 ! Compute all missing vertical scale factor and depths … … 612 613 ! Horizontal scale factor interpolations 613 614 ! -------------------------------------- 614 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt615 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 615 616 ! - JC - hu_b, hv_b, hur_b, hvr_b also 616 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 617 618 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 619 617 620 ! Vertical scale factor interpolations 618 ! ------------------------------------ 619 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 620 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 621 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 622 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 623 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 624 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 625 ! t- and w- points depth 626 ! ---------------------- 627 ! set the isf depth as it is in the initial step 628 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 629 fsdepw_n(:,:,1) = 0.0_wp 630 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 631 621 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 622 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 623 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 624 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' ) 625 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 626 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 627 628 ! t- and w- points depth (set the isf depth as it is in the initial step) 629 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 630 gdepw_n(:,:,1) = 0.0_wp 631 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 632 632 DO jk = 2, jpk 633 633 DO jj = 1,jpj … … 636 636 ! 1 for jk = mikt 637 637 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 638 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)639 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &640 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))641 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)638 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 639 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 640 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 641 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 642 642 END DO 643 643 END DO 644 644 END DO 645 645 646 ! Local depth and Inverse of the local depth of the water column at u- and v- points 647 ! ---------------------------------------------------------------------------------- 648 hu (:,:) = hu_a (:,:) 649 hv (:,:) = hv_a (:,:) 650 651 ! Inverse of the local depth 652 hur(:,:) = hur_a(:,:) 653 hvr(:,:) = hvr_a(:,:) 654 655 ! Local depth of the water column at t- points 656 ! -------------------------------------------- 657 ht(:,:) = 0. 658 DO jk = 1, jpkm1 659 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 646 ! Local depth and Inverse of the local depth of the water 647 ! ------------------------------------------------------- 648 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 649 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 650 ! 651 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 652 DO jk = 2, jpkm1 653 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 660 654 END DO 661 655 662 656 ! Write outputs 663 657 ! ============= 664 CALL iom_put( "e3t" , fse3t_n(:,:,:) )665 CALL iom_put( "e3u" , fse3u_n(:,:,:) )666 CALL iom_put( "e3v" , fse3v_n(:,:,:) )667 CALL iom_put( "e3w" , fse3w_n(:,:,:) )668 CALL iom_put( "tpt_dep" , fsde3w_n(:,:,:) )658 CALL iom_put( "e3t", e3t_n(:,:,:) ) 659 CALL iom_put( "e3u", e3u_n(:,:,:) ) 660 CALL iom_put( "e3v", e3v_n(:,:,:) ) 661 CALL iom_put( "e3w", e3w_n(:,:,:) ) 662 CALL iom_put( "tpt_dep", gde3w_n(:,:,:) ) 669 663 IF( iom_use("e3tdef") ) & 670 CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100* tmask(:,:,:) ) ** 2 )664 CALL iom_put( "e3tdef", ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100. * tmask(:,:,:) ) ** 2 ) 671 665 672 666 ! write restart file 673 667 ! ================== 674 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )675 ! 676 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')677 668 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 669 ! 670 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 671 ! 678 672 END SUBROUTINE dom_vvl_sf_swp 679 673 … … 698 692 !!---------------------------------------------------------------------- 699 693 ! 700 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol')694 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 701 695 ! 702 696 IF(ln_wd) THEN … … 785 779 END SELECT 786 780 ! 787 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol')781 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 788 782 ! 789 783 END SUBROUTINE dom_vvl_interpol … … 816 810 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn ) 817 811 ! 818 id1 = iom_varid( numror, ' fse3t_b', ldstop = .FALSE. )819 id2 = iom_varid( numror, ' fse3t_n', ldstop = .FALSE. )812 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 813 id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 820 814 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 821 815 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) … … 825 819 ! ! --------- ! 826 820 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 827 CALL iom_get( numror, jpdom_autoglo, ' fse3t_b', fse3t_b(:,:,:) )828 CALL iom_get( numror, jpdom_autoglo, ' fse3t_n', fse3t_n(:,:,:) )821 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 822 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 829 823 ! needed to restart if land processor not computed 830 IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'824 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 831 825 WHERE ( tmask(:,:,:) == 0.0_wp ) 832 fse3t_n(:,:,:) = e3t_0(:,:,:)833 fse3t_b(:,:,:) = e3t_0(:,:,:)826 e3t_n(:,:,:) = e3t_0(:,:,:) 827 e3t_b(:,:,:) = e3t_0(:,:,:) 834 828 END WHERE 835 829 IF( neuler == 0 ) THEN 836 fse3t_b(:,:,:) = fse3t_n(:,:,:)830 e3t_b(:,:,:) = e3t_n(:,:,:) 837 831 ENDIF 838 832 ELSE IF( id1 > 0 ) THEN 839 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'840 IF(lwp) write(numout,*) ' fse3t_n set equal to fse3t_b.'833 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 834 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 841 835 IF(lwp) write(numout,*) 'neuler is forced to 0' 842 CALL iom_get( numror, jpdom_autoglo, ' fse3t_b', fse3t_b(:,:,:) )843 fse3t_n(:,:,:) = fse3t_b(:,:,:)836 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 837 e3t_n(:,:,:) = e3t_b(:,:,:) 844 838 neuler = 0 845 839 ELSE IF( id2 > 0 ) THEN 846 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'847 IF(lwp) write(numout,*) ' fse3t_b set equal to fse3t_n.'840 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 841 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 848 842 IF(lwp) write(numout,*) 'neuler is forced to 0' 849 CALL iom_get( numror, jpdom_autoglo, ' fse3t_n', fse3t_n(:,:,:) )850 fse3t_b(:,:,:) = fse3t_n(:,:,:)843 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 844 e3t_b(:,:,:) = e3t_n(:,:,:) 851 845 neuler = 0 852 846 ELSE 853 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'847 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 854 848 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 855 849 IF(lwp) write(numout,*) 'neuler is forced to 0' 856 DO jk =1,jpk857 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &858 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)&859 & + e3t_0(:,:,jk)* (1._wp -tmask(:,:,jk))850 DO jk = 1, jpk 851 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 852 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 853 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 860 854 END DO 861 fse3t_b(:,:,:) = fse3t_n(:,:,:)855 e3t_b(:,:,:) = e3t_n(:,:,:) 862 856 neuler = 0 863 857 ENDIF … … 890 884 ! 891 885 ELSE !* Initialize at "rest" 892 fse3t_b(:,:,:) = e3t_0(:,:,:)893 fse3t_n(:,:,:) = e3t_0(:,:,:)886 e3t_b(:,:,:) = e3t_0(:,:,:) 887 e3t_n(:,:,:) = e3t_0(:,:,:) 894 888 sshn(:,:) = 0.0_wp 895 889 … … 922 916 ! ! all cases ! 923 917 ! ! --------- ! 924 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_b', fse3t_b(:,:,:) )925 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_n', fse3t_n(:,:,:) )918 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) ) 919 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) ) 926 920 ! ! ----------------------- ! 927 921 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 1006 1000 ! 1007 1001 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1008 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )1002 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 1009 1003 ! 1010 1004 IF(lwp) THEN ! Print the choice … … 1020 1014 ! 1021 1015 #if defined key_agrif 1022 IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )1016 IF(.NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' ) 1023 1017 #endif 1024 1018 ! … … 1027 1021 !!====================================================================== 1028 1022 END MODULE domvvl 1029 1030 1031 -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5877 r6141 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 … … 40 40 USE closea ! closed seas 41 41 USE c1d ! 1D vertical configuration 42 ! 42 43 USE in_out_manager ! I/O manager 43 44 USE iom ! I/O library … … 75 76 76 77 !! * Substitutions 77 # include "domzgr_substitute.h90"78 78 # include "vectopt_loop_substitute.h90" 79 79 !!---------------------------------------------------------------------- … … 104 104 INTEGER :: ios 105 105 ! 106 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 106 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 107 107 !!---------------------------------------------------------------------- 108 108 ! … … 122 122 WRITE(numout,*) 'dom_zgr : vertical coordinate' 123 123 WRITE(numout,*) '~~~~~~~' 124 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 125 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 126 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 127 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 128 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 124 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 125 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 126 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 127 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 128 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 129 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 129 130 ENDIF 131 132 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 130 133 131 134 ioptio = 0 ! Check Vertical coordinate options … … 157 160 ! 158 161 IF( nprint == 1 .AND. lwp ) THEN 159 WRITE(numout,*) ' MIN val mbathy ', MINVAL(mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )162 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 160 163 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 161 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )162 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL(e3f_0(:,:,:) ), &163 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL(e3v_0(:,:,:) ), &164 & ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL(e3vw_0(:,:,:)), &165 & ' w ', MINVAL(e3w_0(:,:,:) )164 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) 165 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), & 166 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), & 167 & ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), & 168 & ' w ', MINVAL( e3w_0(:,:,:) ) 166 169 167 170 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 168 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )169 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL(e3f_0(:,:,:) ), &170 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL(e3v_0(:,:,:) ), &171 & ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),&172 & ' w ', MAXVAL(e3w_0(:,:,:) )171 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) ) 172 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), & 173 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), & 174 & ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), & 175 & ' w ', MAXVAL( e3w_0(:,:,:) ) 173 176 ENDIF 174 177 ! … … 381 384 !! - bathy : meter bathymetry (in meters) 382 385 !!---------------------------------------------------------------------- 383 INTEGER :: ji, jj, j l, jk ! dummy loop indices386 INTEGER :: ji, jj, jk ! dummy loop indices 384 387 INTEGER :: inum ! temporary logical unit 385 388 INTEGER :: ierror ! error flag … … 543 546 CALL iom_close( inum ) 544 547 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 548 549 ! set grounded point to 0 550 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 551 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 552 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 553 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 554 END WHERE 545 555 END IF 546 556 ! … … 580 590 ! 581 591 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 582 ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin583 IF ( ln_isfcav ) THEN584 WHERE (bathy == risfdep)585 bathy = 0.0_wp ; risfdep = 0.0_wp586 END WHERE587 END IF588 ! end patch589 592 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 590 593 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth … … 676 679 !! - update bathy : meter bathymetry (in meters) 677 680 !!---------------------------------------------------------------------- 678 !!679 681 INTEGER :: ji, jj, jl ! dummy loop indices 680 682 INTEGER :: icompt, ibtest, ikmax ! temporary integers 681 683 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 682 683 684 !!---------------------------------------------------------------------- 684 685 ! … … 777 778 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 778 779 ENDIF 779 780 IF( lwp .AND. nprint == 1 ) THEN ! control print781 WRITE(numout,*)782 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '783 WRITE(numout,*) ' ------------------'784 CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )785 WRITE(numout,*)786 ENDIF787 780 ! 788 781 CALL wrk_dealloc( jpi, jpj, zbathy ) … … 805 798 !! (min value = 1 over land) 806 799 !!---------------------------------------------------------------------- 807 !!808 800 INTEGER :: ji, jj ! dummy loop indices 809 801 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk … … 837 829 END SUBROUTINE zgr_bot_level 838 830 839 SUBROUTINE zgr_top_level 840 !!---------------------------------------------------------------------- 841 !! *** ROUTINE zgr_bot_level *** 831 832 SUBROUTINE zgr_top_level 833 !!---------------------------------------------------------------------- 834 !! *** ROUTINE zgr_top_level *** 842 835 !! 843 836 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) … … 849 842 !! (min value = 1) 850 843 !!---------------------------------------------------------------------- 851 !!852 844 INTEGER :: ji, jj ! dummy loop indices 853 845 REAL(wp), POINTER, DIMENSION(:,:) :: zmik … … 883 875 END SUBROUTINE zgr_top_level 884 876 877 885 878 SUBROUTINE zgr_zco 886 879 !!---------------------------------------------------------------------- 887 880 !! *** ROUTINE zgr_zco *** 888 881 !! 889 !! ** Purpose : define the z-coordinate system882 !! ** Purpose : define the reference z-coordinate system 890 883 !! 891 884 !! ** Method : set 3D coord. arrays to reference 1D array … … 897 890 ! 898 891 DO jk = 1, jpk 899 gdept_0 900 gdepw_0 901 gde p3w_0(:,:,jk) = gdepw_1d(jk)902 e3t_0 903 e3u_0 904 e3v_0 905 e3f_0 906 e3w_0 907 e3uw_0 908 e3vw_0 892 gdept_0(:,:,jk) = gdept_1d(jk) 893 gdepw_0(:,:,jk) = gdepw_1d(jk) 894 gde3w_0(:,:,jk) = gdepw_1d(jk) 895 e3t_0 (:,:,jk) = e3t_1d (jk) 896 e3u_0 (:,:,jk) = e3t_1d (jk) 897 e3v_0 (:,:,jk) = e3t_1d (jk) 898 e3f_0 (:,:,jk) = e3t_1d (jk) 899 e3w_0 (:,:,jk) = e3w_1d (jk) 900 e3uw_0 (:,:,jk) = e3w_1d (jk) 901 e3vw_0 (:,:,jk) = e3w_1d (jk) 909 902 END DO 910 903 ! … … 919 912 !! 920 913 !! ** Purpose : the depth and vertical scale factor in partial step 921 !! z-coordinate case914 !! reference z-coordinate case 922 915 !! 923 916 !! ** Method : Partial steps : computes the 3D vertical scale factors … … 959 952 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 960 953 !!---------------------------------------------------------------------- 961 !!962 954 INTEGER :: ji, jj, jk ! dummy loop indices 963 955 INTEGER :: ik, it, ikb, ikt ! temporary integers 964 LOGICAL :: ll_print ! Allow control print for debugging965 956 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 966 957 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 967 REAL(wp) :: zmax ! Maximum depth968 958 REAL(wp) :: zdiff ! temporary scalar 969 REAL(wp) :: z refdep! temporary scalar959 REAL(wp) :: zmax ! temporary scalar 970 960 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 971 961 !!--------------------------------------------------------------------- … … 973 963 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 974 964 ! 975 CALL wrk_alloc( jpi, jpj, jpk,zprt )965 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 976 966 ! 977 967 IF(lwp) WRITE(numout,*) … … 979 969 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 980 970 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 981 982 ll_print = .FALSE. ! Local variable for debugging983 984 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth985 WRITE(numout,*)986 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'987 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )988 ENDIF989 990 971 991 972 ! bathymetry in level (from bathy_meter) … … 1006 987 END DO 1007 988 1008 IF ( ln_isfcav ) CALL zgr_isf1009 1010 989 ! Scale factors and depth at T- and W-points 1011 990 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1015 994 e3w_0 (:,:,jk) = e3w_1d (jk) 1016 995 END DO 996 997 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 998 IF ( ln_isfcav ) CALL zgr_isf 999 1000 ! Scale factors and depth at T- and W-points 1001 IF ( .NOT. ln_isfcav ) THEN 1002 DO jj = 1, jpj 1003 DO ji = 1, jpi 1004 ik = mbathy(ji,jj) 1005 IF( ik > 0 ) THEN ! ocean point only 1006 ! max ocean level case 1007 IF( ik == jpkm1 ) THEN 1008 zdepwp = bathy(ji,jj) 1009 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1010 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1011 e3t_0(ji,jj,ik ) = ze3tp 1012 e3t_0(ji,jj,ik+1) = ze3tp 1013 e3w_0(ji,jj,ik ) = ze3wp 1014 e3w_0(ji,jj,ik+1) = ze3tp 1015 gdepw_0(ji,jj,ik+1) = zdepwp 1016 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1017 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1018 ! 1019 ELSE ! standard case 1020 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1021 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1022 ENDIF 1023 !gm Bug? check the gdepw_1d 1024 ! ... on ik 1025 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1026 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1027 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1028 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1029 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1030 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1031 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1032 ! ... on ik+1 1033 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1034 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1035 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1036 ENDIF 1037 ENDIF 1038 END DO 1039 END DO 1040 ! 1041 it = 0 1042 DO jj = 1, jpj 1043 DO ji = 1, jpi 1044 ik = mbathy(ji,jj) 1045 IF( ik > 0 ) THEN ! ocean point only 1046 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1047 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1048 ! test 1049 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1050 IF( zdiff <= 0._wp .AND. lwp ) THEN 1051 it = it + 1 1052 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1053 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1054 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1055 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1056 ENDIF 1057 ENDIF 1058 END DO 1059 END DO 1060 END IF 1061 ! 1062 ! Scale factors and depth at U-, V-, UW and VW-points 1063 DO jk = 1, jpk ! initialisation to z-scale factors 1064 e3u_0 (:,:,jk) = e3t_1d(jk) 1065 e3v_0 (:,:,jk) = e3t_1d(jk) 1066 e3uw_0(:,:,jk) = e3w_1d(jk) 1067 e3vw_0(:,:,jk) = e3w_1d(jk) 1068 END DO 1069 1070 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1071 DO jj = 1, jpjm1 1072 DO ji = 1, fs_jpim1 ! vector opt. 1073 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1074 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1075 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1076 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1077 END DO 1078 END DO 1079 END DO 1080 IF ( ln_isfcav ) THEN 1081 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1082 DO jj = 2, jpjm1 1083 DO ji = 2, fs_jpim1 ! vector opt. 1084 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1085 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1086 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1087 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1088 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1089 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1090 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1091 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1092 END DO 1093 END DO 1094 END IF 1095 1096 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1097 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1098 ! 1099 1100 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1101 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1102 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1103 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1104 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1105 END DO 1106 1107 ! Scale factor at F-point 1108 DO jk = 1, jpk ! initialisation to z-scale factors 1109 e3f_0(:,:,jk) = e3t_1d(jk) 1110 END DO 1111 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1112 DO jj = 1, jpjm1 1113 DO ji = 1, fs_jpim1 ! vector opt. 1114 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1115 END DO 1116 END DO 1117 END DO 1118 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1119 ! 1120 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1121 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1122 END DO 1123 !!gm bug ? : must be a do loop with mj0,mj1 1017 1124 ! 1125 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1126 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1127 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1128 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1129 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1130 1131 ! Control of the sign 1132 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1133 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1134 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1135 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1136 1137 ! Compute gde3w_0 (vertical sum of e3w) 1138 IF ( ln_isfcav ) THEN ! if cavity 1139 WHERE( misfdep == 0 ) misfdep = 1 1140 DO jj = 1,jpj 1141 DO ji = 1,jpi 1142 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1143 DO jk = 2, misfdep(ji,jj) 1144 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1145 END DO 1146 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)) 1147 DO jk = misfdep(ji,jj) + 1, jpk 1148 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1149 END DO 1150 END DO 1151 END DO 1152 ELSE ! no cavity 1153 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1154 DO jk = 2, jpk 1155 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1156 END DO 1157 END IF 1158 ! 1159 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1160 ! 1161 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1162 ! 1163 END SUBROUTINE zgr_zps 1164 1165 1166 SUBROUTINE zgr_isf 1167 !!---------------------------------------------------------------------- 1168 !! *** ROUTINE zgr_isf *** 1169 !! 1170 !! ** Purpose : check the bathymetry in levels 1171 !! 1172 !! ** Method : THe water column have to contained at least 2 cells 1173 !! Bathymetry and isfdraft are modified (dig/close) to respect 1174 !! this criterion. 1175 !! 1176 !! ** Action : - test compatibility between isfdraft and bathy 1177 !! - bathy and isfdraft are modified 1178 !!---------------------------------------------------------------------- 1179 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1180 INTEGER :: ik, it ! temporary integers 1181 INTEGER :: icompt, ibtest ! (ISF) 1182 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1183 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1184 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1185 REAL(wp) :: zmax ! Maximum and minimum depth 1186 REAL(wp) :: zbathydiff ! isf temporary scalar 1187 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1188 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1189 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1190 REAL(wp) :: zdiff ! temporary scalar 1191 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1192 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1193 !!--------------------------------------------------------------------- 1194 ! 1195 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1196 ! 1197 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1198 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1199 1200 ! (ISF) compute misfdep 1201 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1202 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1203 END WHERE 1204 1205 ! Compute misfdep for ocean points (i.e. first wet level) 1206 ! find the first ocean level such that the first level thickness 1207 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1208 ! e3t_0 is the reference level thickness 1209 DO jk = 2, jpkm1 1210 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1211 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1212 END DO 1213 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1214 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1215 END WHERE 1216 1217 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1218 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1219 misfdep = 0; risfdep = 0.0_wp; 1220 mbathy = 0; bathy = 0.0_wp; 1221 END WHERE 1222 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1223 misfdep = 0; risfdep = 0.0_wp; 1224 mbathy = 0; bathy = 0.0_wp; 1225 END WHERE 1226 1227 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 1228 icompt = 0 1229 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1230 DO jl = 1, 10 1231 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 1232 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 1233 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1234 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1235 END WHERE 1236 WHERE (mbathy(:,:) <= 0) 1237 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1238 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1239 END WHERE 1240 IF( lk_mpp ) THEN 1241 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1242 CALL lbc_lnk( zbathy, 'T', 1. ) 1243 misfdep(:,:) = INT( zbathy(:,:) ) 1244 1245 CALL lbc_lnk( risfdep,'T', 1. ) 1246 CALL lbc_lnk( bathy, 'T', 1. ) 1247 1248 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1249 CALL lbc_lnk( zbathy, 'T', 1. ) 1250 mbathy(:,:) = INT( zbathy(:,:) ) 1251 ENDIF 1252 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1253 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1254 misfdep(jpi,:) = misfdep( 2 ,:) 1255 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1256 mbathy(jpi,:) = mbathy( 2 ,:) 1257 ENDIF 1258 1259 ! split last cell if possible (only where water column is 2 cell or less) 1260 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 1261 IF ( .NOT. ln_iscpl) THEN 1262 DO jk = jpkm1, 1, -1 1263 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1264 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1265 mbathy(:,:) = jk 1266 bathy(:,:) = zmax 1267 END WHERE 1268 END DO 1269 END IF 1270 1271 ! split top cell if possible (only where water column is 2 cell or less) 1272 DO jk = 2, jpkm1 1273 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1274 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1275 misfdep(:,:) = jk 1276 risfdep(:,:) = zmax 1277 END WHERE 1278 END DO 1279 1280 1281 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1282 DO jj = 1, jpj 1283 DO ji = 1, jpi 1284 ! find the minimum change option: 1285 ! test bathy 1286 IF (risfdep(ji,jj) > 1) THEN 1287 IF ( .NOT. ln_iscpl ) THEN 1288 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1289 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1290 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1291 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1292 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1293 IF (zbathydiff <= zrisfdepdiff) THEN 1294 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1295 mbathy(ji,jj)= mbathy(ji,jj) + 1 1296 ELSE 1297 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1298 misfdep(ji,jj) = misfdep(ji,jj) - 1 1299 END IF 1300 ENDIF 1301 ELSE 1302 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1303 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1304 misfdep(ji,jj) = misfdep(ji,jj) - 1 1305 END IF 1306 END IF 1307 END IF 1308 END DO 1309 END DO 1310 1311 ! At least 2 levels for water thickness at T, U, and V point. 1312 DO jj = 1, jpj 1313 DO ji = 1, jpi 1314 ! find the minimum change option: 1315 ! test bathy 1316 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1317 IF ( .NOT. ln_iscpl ) THEN 1318 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1319 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1320 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) & 1321 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1322 IF (zbathydiff <= zrisfdepdiff) THEN 1323 mbathy(ji,jj) = mbathy(ji,jj) + 1 1324 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1325 ELSE 1326 misfdep(ji,jj)= misfdep(ji,jj) - 1 1327 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1328 END IF 1329 ELSE 1330 misfdep(ji,jj)= misfdep(ji,jj) - 1 1331 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1332 END IF 1333 ENDIF 1334 END DO 1335 END DO 1336 1337 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1338 DO jj = 1, jpjm1 1339 DO ji = 1, jpim1 1340 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1341 IF ( .NOT. ln_iscpl ) THEN 1342 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) & 1343 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1344 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 1345 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1346 IF (zbathydiff <= zrisfdepdiff) THEN 1347 mbathy(ji,jj) = mbathy(ji,jj) + 1 1348 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1349 ELSE 1350 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1351 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1352 END IF 1353 ELSE 1354 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1355 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1356 END IF 1357 ENDIF 1358 END DO 1359 END DO 1360 1361 IF( lk_mpp ) THEN 1362 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1363 CALL lbc_lnk( zbathy, 'T', 1. ) 1364 misfdep(:,:) = INT( zbathy(:,:) ) 1365 1366 CALL lbc_lnk( risfdep,'T', 1. ) 1367 CALL lbc_lnk( bathy, 'T', 1. ) 1368 1369 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1370 CALL lbc_lnk( zbathy, 'T', 1. ) 1371 mbathy(:,:) = INT( zbathy(:,:) ) 1372 ENDIF 1373 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1374 DO jj = 1, jpjm1 1375 DO ji = 1, jpim1 1376 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1377 IF ( .NOT. ln_iscpl ) THEN 1378 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 1379 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1380 zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) & 1381 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1382 IF (zbathydiff <= zrisfdepdiff) THEN 1383 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1384 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1385 ELSE 1386 misfdep(ji,jj) = misfdep(ji,jj) - 1 1387 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1388 END IF 1389 ELSE 1390 misfdep(ji,jj) = misfdep(ji,jj) - 1 1391 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1392 END IF 1393 ENDIF 1394 END DO 1395 END DO 1396 1397 1398 IF( lk_mpp ) THEN 1399 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1400 CALL lbc_lnk( zbathy, 'T', 1. ) 1401 misfdep(:,:) = INT( zbathy(:,:) ) 1402 1403 CALL lbc_lnk( risfdep,'T', 1. ) 1404 CALL lbc_lnk( bathy, 'T', 1. ) 1405 1406 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1407 CALL lbc_lnk( zbathy, 'T', 1. ) 1408 mbathy(:,:) = INT( zbathy(:,:) ) 1409 ENDIF 1410 1411 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1412 DO jj = 1, jpjm1 1413 DO ji = 1, jpim1 1414 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1415 IF ( .NOT. ln_iscpl ) THEN 1416 zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1417 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1418 zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 1419 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1420 IF (zbathydiff <= zrisfdepdiff) THEN 1421 mbathy(ji,jj) = mbathy(ji,jj) + 1 1422 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1423 ELSE 1424 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1425 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1426 END IF 1427 ELSE 1428 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1429 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1430 ENDIF 1431 ENDIF 1432 ENDDO 1433 ENDDO 1434 1435 IF( lk_mpp ) THEN 1436 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1437 CALL lbc_lnk( zbathy, 'T', 1. ) 1438 misfdep(:,:) = INT( zbathy(:,:) ) 1439 1440 CALL lbc_lnk( risfdep,'T', 1. ) 1441 CALL lbc_lnk( bathy, 'T', 1. ) 1442 1443 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1444 CALL lbc_lnk( zbathy, 'T', 1. ) 1445 mbathy(:,:) = INT( zbathy(:,:) ) 1446 ENDIF 1447 1448 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1449 DO jj = 1, jpjm1 1450 DO ji = 1, jpim1 1451 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1452 IF ( .NOT. ln_iscpl ) THEN 1453 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 1454 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1455 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) & 1456 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1457 IF (zbathydiff <= zrisfdepdiff) THEN 1458 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1459 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1460 ELSE 1461 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1462 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1463 END IF 1464 ELSE 1465 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1466 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1467 ENDIF 1468 ENDIF 1469 ENDDO 1470 ENDDO 1471 1472 IF( lk_mpp ) THEN 1473 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1474 CALL lbc_lnk( zbathy, 'T', 1. ) 1475 misfdep(:,:) = INT( zbathy(:,:) ) 1476 1477 CALL lbc_lnk( risfdep,'T', 1. ) 1478 CALL lbc_lnk( bathy, 'T', 1. ) 1479 1480 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1481 CALL lbc_lnk( zbathy, 'T', 1. ) 1482 mbathy(:,:) = INT( zbathy(:,:) ) 1483 ENDIF 1484 END DO 1485 ! end dig bathy/ice shelf to be compatible 1486 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1487 DO jl = 1,20 1488 1489 ! remove single point "bay" on isf coast line in the ice shelf draft' 1490 DO jk = 2, jpk 1491 WHERE (misfdep==0) misfdep=jpk 1492 zmask=0._wp 1493 WHERE (misfdep <= jk) zmask=1 1494 DO jj = 2, jpjm1 1495 DO ji = 2, jpim1 1496 IF (misfdep(ji,jj) == jk) THEN 1497 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1498 IF (ibtest <= 1) THEN 1499 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1500 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1501 END IF 1502 END IF 1503 END DO 1504 END DO 1505 END DO 1506 WHERE (misfdep==jpk) 1507 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1508 END WHERE 1509 IF( lk_mpp ) THEN 1510 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1511 CALL lbc_lnk( zbathy, 'T', 1. ) 1512 misfdep(:,:) = INT( zbathy(:,:) ) 1513 1514 CALL lbc_lnk( risfdep,'T', 1. ) 1515 CALL lbc_lnk( bathy, 'T', 1. ) 1516 1517 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1518 CALL lbc_lnk( zbathy, 'T', 1. ) 1519 mbathy(:,:) = INT( zbathy(:,:) ) 1520 ENDIF 1521 1522 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1523 DO jk = jpk,1,-1 1524 zmask=0._wp 1525 WHERE (mbathy >= jk ) zmask=1 1526 DO jj = 2, jpjm1 1527 DO ji = 2, jpim1 1528 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1529 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1530 IF (ibtest <= 1) THEN 1531 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1532 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1533 END IF 1534 END IF 1535 END DO 1536 END DO 1537 END DO 1538 WHERE (mbathy==0) 1539 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1540 END WHERE 1541 IF( lk_mpp ) THEN 1542 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1543 CALL lbc_lnk( zbathy, 'T', 1. ) 1544 misfdep(:,:) = INT( zbathy(:,:) ) 1545 1546 CALL lbc_lnk( risfdep,'T', 1. ) 1547 CALL lbc_lnk( bathy, 'T', 1. ) 1548 1549 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1550 CALL lbc_lnk( zbathy, 'T', 1. ) 1551 mbathy(:,:) = INT( zbathy(:,:) ) 1552 ENDIF 1553 1554 ! fill hole in ice shelf 1555 zmisfdep = misfdep 1556 zrisfdep = risfdep 1557 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1558 DO jj = 2, jpjm1 1559 DO ji = 2, jpim1 1560 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1561 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1562 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk 1563 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk 1564 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk 1565 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk 1566 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1567 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1568 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1569 END IF 1570 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1571 misfdep(ji,jj) = ibtest 1572 risfdep(ji,jj) = gdepw_1d(ibtest) 1573 ENDIF 1574 ENDDO 1575 ENDDO 1576 1577 IF( lk_mpp ) THEN 1578 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1579 CALL lbc_lnk( zbathy, 'T', 1. ) 1580 misfdep(:,:) = INT( zbathy(:,:) ) 1581 1582 CALL lbc_lnk( risfdep, 'T', 1. ) 1583 CALL lbc_lnk( bathy, 'T', 1. ) 1584 1585 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1586 CALL lbc_lnk( zbathy, 'T', 1. ) 1587 mbathy(:,:) = INT( zbathy(:,:) ) 1588 ENDIF 1589 ! 1590 !! fill hole in bathymetry 1591 zmbathy (:,:)=mbathy (:,:) 1592 DO jj = 2, jpjm1 1593 DO ji = 2, jpim1 1594 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1595 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1596 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0 1597 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1598 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1599 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1600 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1601 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1602 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1603 END IF 1604 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1605 mbathy(ji,jj) = ibtest 1606 bathy(ji,jj) = gdepw_1d(ibtest+1) 1607 ENDIF 1608 END DO 1609 END DO 1610 IF( lk_mpp ) THEN 1611 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1612 CALL lbc_lnk( zbathy, 'T', 1. ) 1613 misfdep(:,:) = INT( zbathy(:,:) ) 1614 1615 CALL lbc_lnk( risfdep, 'T', 1. ) 1616 CALL lbc_lnk( bathy, 'T', 1. ) 1617 1618 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1619 CALL lbc_lnk( zbathy, 'T', 1. ) 1620 mbathy(:,:) = INT( zbathy(:,:) ) 1621 ENDIF 1622 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1623 DO jj = 1, jpjm1 1624 DO ji = 1, jpim1 1625 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1626 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1627 END IF 1628 END DO 1629 END DO 1630 IF( lk_mpp ) THEN 1631 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1632 CALL lbc_lnk( zbathy, 'T', 1. ) 1633 misfdep(:,:) = INT( zbathy(:,:) ) 1634 1635 CALL lbc_lnk( risfdep, 'T', 1. ) 1636 CALL lbc_lnk( bathy, 'T', 1. ) 1637 1638 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1639 CALL lbc_lnk( zbathy, 'T', 1. ) 1640 mbathy(:,:) = INT( zbathy(:,:) ) 1641 ENDIF 1642 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1643 DO jj = 1, jpjm1 1644 DO ji = 1, jpim1 1645 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1646 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1647 END IF 1648 END DO 1649 END DO 1650 IF( lk_mpp ) THEN 1651 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1652 CALL lbc_lnk( zbathy, 'T', 1. ) 1653 misfdep(:,:) = INT( zbathy(:,:) ) 1654 1655 CALL lbc_lnk( risfdep,'T', 1. ) 1656 CALL lbc_lnk( bathy, 'T', 1. ) 1657 1658 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1659 CALL lbc_lnk( zbathy, 'T', 1. ) 1660 mbathy(:,:) = INT( zbathy(:,:) ) 1661 ENDIF 1662 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1663 DO jj = 1, jpjm1 1664 DO ji = 1, jpi 1665 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1666 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1667 END IF 1668 END DO 1669 END DO 1670 IF( lk_mpp ) THEN 1671 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1672 CALL lbc_lnk( zbathy, 'T', 1. ) 1673 misfdep(:,:) = INT( zbathy(:,:) ) 1674 1675 CALL lbc_lnk( risfdep,'T', 1. ) 1676 CALL lbc_lnk( bathy, 'T', 1. ) 1677 1678 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1679 CALL lbc_lnk( zbathy, 'T', 1. ) 1680 mbathy(:,:) = INT( zbathy(:,:) ) 1681 ENDIF 1682 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1683 DO jj = 1, jpjm1 1684 DO ji = 1, jpi 1685 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1686 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1687 END IF 1688 END DO 1689 END DO 1690 IF( lk_mpp ) THEN 1691 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1692 CALL lbc_lnk( zbathy, 'T', 1. ) 1693 misfdep(:,:) = INT( zbathy(:,:) ) 1694 1695 CALL lbc_lnk( risfdep,'T', 1. ) 1696 CALL lbc_lnk( bathy, 'T', 1. ) 1697 1698 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1699 CALL lbc_lnk( zbathy, 'T', 1. ) 1700 mbathy(:,:) = INT( zbathy(:,:) ) 1701 ENDIF 1702 ! if not compatible after all check, mask T 1703 DO jj = 1, jpj 1704 DO ji = 1, jpi 1705 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1706 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1707 END IF 1708 END DO 1709 END DO 1710 1711 WHERE (mbathy(:,:) == 1) 1712 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1713 END WHERE 1714 END DO 1715 ! end check compatibility ice shelf/bathy 1716 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1717 WHERE (risfdep(:,:) <= 10._wp) 1718 misfdep = 1; risfdep = 0.0_wp; 1719 END WHERE 1720 1721 IF( icompt == 0 ) THEN 1722 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1723 ELSE 1724 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1725 ENDIF 1726 1727 ! compute scale factor and depth at T- and W- points 1018 1728 DO jj = 1, jpj 1019 1729 DO ji = 1, jpi … … 1037 1747 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1038 1748 ENDIF 1749 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1039 1750 !gm Bug? check the gdepw_1d 1040 1751 ! ... on ik … … 1042 1753 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1043 1754 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1044 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1045 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1046 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1047 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1755 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1756 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1048 1757 ! ... on ik+1 1049 1758 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1050 1759 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1051 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1052 1760 ENDIF 1053 1761 ENDIF … … 1075 1783 END DO 1076 1784 ! 1077 IF ( ln_isfcav ) THEN1078 1785 ! (ISF) Definition of e3t, u, v, w for ISF case 1079 1080 1081 1082 1083 1084 1786 DO jj = 1, jpj 1787 DO ji = 1, jpi 1788 ik = misfdep(ji,jj) 1789 IF( ik > 1 ) THEN ! ice shelf point only 1790 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1791 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1085 1792 !gm Bug? check the gdepw_0 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 e3w_0 (ji,jj,ik ) =2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))1098 1793 ! ... on ik 1794 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1795 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1796 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1797 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1798 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1799 1800 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1801 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1802 ENDIF 1803 ! ... on ik / ik-1 1804 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1805 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1099 1806 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1100 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1807 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1808 ENDIF 1809 END DO 1810 END DO 1811 1812 it = 0 1813 DO jj = 1, jpj 1814 DO ji = 1, jpi 1815 ik = misfdep(ji,jj) 1816 IF( ik > 1 ) THEN ! ice shelf point only 1817 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1818 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1819 ! test 1820 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1821 IF( zdiff <= 0. .AND. lwp ) THEN 1822 it = it + 1 1823 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1824 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1825 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1826 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1101 1827 ENDIF 1102 END DO1828 ENDIF 1103 1829 END DO 1104 !1105 it = 01106 DO jj = 1, jpj1107 DO ji = 1, jpi1108 ik = misfdep(ji,jj)1109 IF( ik > 1 ) THEN ! ice shelf point only1110 e3tp (ji,jj) = e3t_0(ji,jj,ik )1111 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )1112 ! test1113 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )1114 IF( zdiff <= 0. .AND. lwp ) THEN1115 it = it + 11116 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1117 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)1118 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1119 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)1120 ENDIF1121 ENDIF1122 END DO1123 END DO1124 END IF1125 ! END (ISF)1126 1127 ! Scale factors and depth at U-, V-, UW and VW-points1128 DO jk = 1, jpk ! initialisation to z-scale factors1129 e3u_0 (:,:,jk) = e3t_1d(jk)1130 e3v_0 (:,:,jk) = e3t_1d(jk)1131 e3uw_0(:,:,jk) = e3w_1d(jk)1132 e3vw_0(:,:,jk) = e3w_1d(jk)1133 END DO1134 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors1135 DO jj = 1, jpjm11136 DO ji = 1, fs_jpim1 ! vector opt.1137 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )1138 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )1139 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )1140 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )1141 END DO1142 END DO1143 END DO1144 IF ( ln_isfcav ) THEN1145 ! (ISF) define e3uw (adapted for 2 cells in the water column)1146 DO jj = 2, jpjm11147 DO ji = 2, fs_jpim1 ! vector opt.1148 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))1149 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))1150 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &1151 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )1152 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))1153 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))1154 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &1155 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )1156 END DO1157 END DO1158 END IF1159 1160 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions1161 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )1162 !1163 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1164 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)1165 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)1166 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)1167 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)1168 END DO1169 1170 ! Scale factor at F-point1171 DO jk = 1, jpk ! initialisation to z-scale factors1172 e3f_0(:,:,jk) = e3t_1d(jk)1173 END DO1174 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors1175 DO jj = 1, jpjm11176 DO ji = 1, fs_jpim1 ! vector opt.1177 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )1178 END DO1179 END DO1180 END DO1181 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions1182 !1183 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)1184 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)1185 END DO1186 !!gm bug ? : must be a do loop with mj0,mj11187 !1188 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 21189 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)1190 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)1191 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)1192 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)1193 1194 ! Control of the sign1195 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )1196 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )1197 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )1198 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )1199 1200 ! Compute gdep3w_0 (vertical sum of e3w)1201 IF ( ln_isfcav ) THEN ! if cavity1202 WHERE (misfdep == 0) misfdep = 11203 DO jj = 1,jpj1204 DO ji = 1,jpi1205 gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1206 DO jk = 2, misfdep(ji,jj)1207 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1208 END DO1209 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))1210 DO jk = misfdep(ji,jj) + 1, jpk1211 gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1212 END DO1213 END DO1214 END DO1215 ELSE ! no cavity1216 gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1217 DO jk = 2, jpk1218 gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)1219 END DO1220 END IF1221 ! ! ================= !1222 IF(lwp .AND. ll_print) THEN ! Control print !1223 ! ! ================= !1224 DO jj = 1,jpj1225 DO ji = 1, jpi1226 ik = MAX( mbathy(ji,jj), 1 )1227 zprt(ji,jj,1) = e3t_0 (ji,jj,ik)1228 zprt(ji,jj,2) = e3w_0 (ji,jj,ik)1229 zprt(ji,jj,3) = e3u_0 (ji,jj,ik)1230 zprt(ji,jj,4) = e3v_0 (ji,jj,ik)1231 zprt(ji,jj,5) = e3f_0 (ji,jj,ik)1232 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)1233 END DO1234 END DO1235 WRITE(numout,*)1236 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1237 WRITE(numout,*)1238 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1239 WRITE(numout,*)1240 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1241 WRITE(numout,*)1242 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1243 WRITE(numout,*)1244 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1245 WRITE(numout,*)1246 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)1247 ENDIF1248 !1249 CALL wrk_dealloc( jpi, jpj, jpk, zprt )1250 !1251 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')1252 !1253 END SUBROUTINE zgr_zps1254 1255 SUBROUTINE zgr_isf1256 !!----------------------------------------------------------------------1257 !! *** ROUTINE zgr_isf ***1258 !!1259 !! ** Purpose : check the bathymetry in levels1260 !!1261 !! ** Method : THe water column have to contained at least 2 cells1262 !! Bathymetry and isfdraft are modified (dig/close) to respect1263 !! this criterion.1264 !!1265 !!1266 !! ** Action : - test compatibility between isfdraft and bathy1267 !! - bathy and isfdraft are modified1268 !!----------------------------------------------------------------------1269 !!1270 INTEGER :: ji, jj, jk, jl ! dummy loop indices1271 INTEGER :: ik, it ! temporary integers1272 INTEGER :: id, jd, nprocd1273 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)1274 LOGICAL :: ll_print ! Allow control print for debugging1275 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points1276 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t1277 REAL(wp) :: zmax, zmin ! Maximum and minimum depth1278 REAL(wp) :: zdiff ! temporary scalar1279 REAL(wp) :: zrefdep ! temporary scalar1280 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1281 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)1282 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)1283 !!---------------------------------------------------------------------1284 !1285 IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1286 !1287 CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)1288 CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )1289 1290 1291 ! (ISF) compute misfdep1292 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 11293 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level1294 END WHERE1295 1296 ! Compute misfdep for ocean points (i.e. first wet level)1297 ! find the first ocean level such that the first level thickness1298 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where1299 ! e3t_0 is the reference level thickness1300 DO jk = 2, jpkm11301 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1302 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+11303 1830 END DO 1304 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)1305 risfdep(:,:) = 0. ; misfdep(:,:) = 11306 END WHERE1307 1308 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation1309 icompt = 01310 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together1311 DO jl = 1, 101312 WHERE (bathy(:,:) .EQ. risfdep(:,:) )1313 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp1314 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp1315 END WHERE1316 WHERE (mbathy(:,:) <= 0)1317 misfdep(:,:) = 0; risfdep(:,:) = 0._wp1318 mbathy (:,:) = 0; bathy (:,:) = 0._wp1319 ENDWHERE1320 IF( lk_mpp ) THEN1321 zbathy(:,:) = FLOAT( misfdep(:,:) )1322 CALL lbc_lnk( zbathy, 'T', 1. )1323 misfdep(:,:) = INT( zbathy(:,:) )1324 CALL lbc_lnk( risfdep, 'T', 1. )1325 CALL lbc_lnk( bathy, 'T', 1. )1326 zbathy(:,:) = FLOAT( mbathy(:,:) )1327 CALL lbc_lnk( zbathy, 'T', 1. )1328 mbathy(:,:) = INT( zbathy(:,:) )1329 ENDIF1330 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1331 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west1332 misfdep(jpi,:) = misfdep( 2 ,:)1333 ENDIF1334 1335 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN1336 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west1337 mbathy(jpi,:) = mbathy( 2 ,:)1338 ENDIF1339 1340 ! split last cell if possible (only where water column is 2 cell or less)1341 DO jk = jpkm1, 1, -11342 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1343 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)1344 mbathy(:,:) = jk1345 bathy(:,:) = zmax1346 END WHERE1347 END DO1348 1349 ! split top cell if possible (only where water column is 2 cell or less)1350 DO jk = 2, jpkm11351 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )1352 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)1353 misfdep(:,:) = jk1354 risfdep(:,:) = zmax1355 END WHERE1356 END DO1357 1358 1359 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition1360 DO jj = 1, jpj1361 DO ji = 1, jpi1362 ! find the minimum change option:1363 ! test bathy1364 IF (risfdep(ji,jj) .GT. 1) THEN1365 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1366 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1367 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1368 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1369 1370 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN1371 IF (zbathydiff .LE. zrisfdepdiff) THEN1372 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )1373 mbathy(ji,jj)= mbathy(ji,jj) + 11374 ELSE1375 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )1376 misfdep(ji,jj) = misfdep(ji,jj) - 11377 END IF1378 END IF1379 END IF1380 END DO1381 END DO1382 1383 ! At least 2 levels for water thickness at T, U, and V point.1384 DO jj = 1, jpj1385 DO ji = 1, jpi1386 ! find the minimum change option:1387 ! test bathy1388 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1389 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1)&1390 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))1391 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &1392 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))1393 IF (zbathydiff .LE. zrisfdepdiff) THEN1394 mbathy(ji,jj) = mbathy(ji,jj) + 11395 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1396 ELSE1397 misfdep(ji,jj)= misfdep(ji,jj) - 11398 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )1399 END IF1400 ENDIF1401 END DO1402 END DO1403 1404 ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)1405 DO jj = 1, jpjm11406 DO ji = 1, jpim11407 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1408 zbathydiff =ABS(bathy(ji,jj ) - (gdepw_1d(mbathy (ji,jj)+1) &1409 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))1410 zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) &1411 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))1412 IF (zbathydiff .LE. zrisfdepdiff) THEN1413 mbathy(ji,jj) = mbathy(ji,jj) + 11414 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) &1415 & + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )1416 ELSE1417 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 11418 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &1419 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )1420 END IF1421 ENDIF1422 END DO1423 END DO1424 1425 IF( lk_mpp ) THEN1426 zbathy(:,:) = FLOAT( misfdep(:,:) )1427 CALL lbc_lnk( zbathy, 'T', 1. )1428 misfdep(:,:) = INT( zbathy(:,:) )1429 CALL lbc_lnk( risfdep, 'T', 1. )1430 CALL lbc_lnk( bathy, 'T', 1. )1431 zbathy(:,:) = FLOAT( mbathy(:,:) )1432 CALL lbc_lnk( zbathy, 'T', 1. )1433 mbathy(:,:) = INT( zbathy(:,:) )1434 ENDIF1435 ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)1436 DO jj = 1, jpjm11437 DO ji = 1, jpim11438 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN1439 zbathydiff =ABS( bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &1440 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))1441 zrisfdepdiff=ABS(risfdep(ji,jj ) - (gdepw_1d(misfdep(ji,jj ) ) &1442 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))1443 IF (zbathydiff .LE. zrisfdepdiff) THEN1444 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 11445 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )1446 ELSE1447 misfdep(ji,jj) = misfdep(ji,jj) - 11448 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )1449 END IF1450 ENDIF1451 END DO1452 END DO1453 1454 1455 IF( lk_mpp ) THEN1456 zbathy(:,:) = FLOAT( misfdep(:,:) )1457 CALL lbc_lnk( zbathy, 'T', 1. )1458 misfdep(:,:) = INT( zbathy(:,:) )1459 CALL lbc_lnk( risfdep, 'T', 1. )1460 CALL lbc_lnk( bathy, 'T', 1. )1461 zbathy(:,:) = FLOAT( mbathy(:,:) )1462 CALL lbc_lnk( zbathy, 'T', 1. )1463 mbathy(:,:) = INT( zbathy(:,:) )1464 ENDIF1465 1466 ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)1467 DO jj = 1, jpjm11468 DO ji = 1, jpim11469 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1470 zbathydiff =ABS( bathy(ji ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &1471 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))1472 zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &1473 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))1474 IF (zbathydiff .LE. zrisfdepdiff) THEN1475 mbathy(ji,jj) = mbathy(ji,jj) + 11476 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )1477 ELSE1478 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 11479 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )1480 END IF1481 ENDIF1482 ENDDO1483 ENDDO1484 1485 IF( lk_mpp ) THEN1486 zbathy(:,:) = FLOAT( misfdep(:,:) )1487 CALL lbc_lnk( zbathy, 'T', 1. )1488 misfdep(:,:) = INT( zbathy(:,:) )1489 CALL lbc_lnk( risfdep, 'T', 1. )1490 CALL lbc_lnk( bathy, 'T', 1. )1491 zbathy(:,:) = FLOAT( mbathy(:,:) )1492 CALL lbc_lnk( zbathy, 'T', 1. )1493 mbathy(:,:) = INT( zbathy(:,:) )1494 ENDIF1495 1496 ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)1497 DO jj = 1, jpjm11498 DO ji = 1, jpim11499 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN1500 zbathydiff =ABS( bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &1501 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))1502 zrisfdepdiff=ABS(risfdep(ji ,jj) - (gdepw_1d(misfdep(ji ,jj) ) &1503 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))1504 IF (zbathydiff .LE. zrisfdepdiff) THEN1505 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 11506 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) &1507 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )1508 ELSE1509 misfdep(ji,jj) = misfdep(ji ,jj) - 11510 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) &1511 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )1512 END IF1513 ENDIF1514 ENDDO1515 ENDDO1516 1517 IF( lk_mpp ) THEN1518 zbathy(:,:) = FLOAT( misfdep(:,:) )1519 CALL lbc_lnk( zbathy, 'T', 1. )1520 misfdep(:,:) = INT( zbathy(:,:) )1521 CALL lbc_lnk( risfdep, 'T', 1. )1522 CALL lbc_lnk( bathy, 'T', 1. )1523 zbathy(:,:) = FLOAT( mbathy(:,:) )1524 CALL lbc_lnk( zbathy, 'T', 1. )1525 mbathy(:,:) = INT( zbathy(:,:) )1526 ENDIF1527 END DO1528 ! end dig bathy/ice shelf to be compatible1529 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness1530 DO jl = 1,201531 1532 ! remove single point "bay" on isf coast line in the ice shelf draft'1533 DO jk = 2, jpk1534 WHERE (misfdep==0) misfdep=jpk1535 zmask=01536 WHERE (misfdep .LE. jk) zmask=11537 DO jj = 2, jpjm11538 DO ji = 2, jpim11539 IF (misfdep(ji,jj) .EQ. jk) THEN1540 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1541 IF (ibtest .LE. 1) THEN1542 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+11543 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk1544 END IF1545 END IF1546 END DO1547 END DO1548 END DO1549 WHERE (misfdep==jpk)1550 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1551 END WHERE1552 IF( lk_mpp ) THEN1553 zbathy(:,:) = FLOAT( misfdep(:,:) )1554 CALL lbc_lnk( zbathy, 'T', 1. )1555 misfdep(:,:) = INT( zbathy(:,:) )1556 CALL lbc_lnk( risfdep, 'T', 1. )1557 CALL lbc_lnk( bathy, 'T', 1. )1558 zbathy(:,:) = FLOAT( mbathy(:,:) )1559 CALL lbc_lnk( zbathy, 'T', 1. )1560 mbathy(:,:) = INT( zbathy(:,:) )1561 ENDIF1562 1563 ! remove single point "bay" on bathy coast line beneath an ice shelf'1564 DO jk = jpk,1,-11565 zmask=01566 WHERE (mbathy .GE. jk ) zmask=11567 DO jj = 2, jpjm11568 DO ji = 2, jpim11569 IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN1570 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)1571 IF (ibtest .LE. 1) THEN1572 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-11573 IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 01574 END IF1575 END IF1576 END DO1577 END DO1578 END DO1579 WHERE (mbathy==0)1580 misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.1581 END WHERE1582 IF( lk_mpp ) THEN1583 zbathy(:,:) = FLOAT( misfdep(:,:) )1584 CALL lbc_lnk( zbathy, 'T', 1. )1585 misfdep(:,:) = INT( zbathy(:,:) )1586 CALL lbc_lnk( risfdep, 'T', 1. )1587 CALL lbc_lnk( bathy, 'T', 1. )1588 zbathy(:,:) = FLOAT( mbathy(:,:) )1589 CALL lbc_lnk( zbathy, 'T', 1. )1590 mbathy(:,:) = INT( zbathy(:,:) )1591 ENDIF1592 1593 ! fill hole in ice shelf1594 zmisfdep = misfdep1595 zrisfdep = risfdep1596 WHERE (zmisfdep .LE. 1) zmisfdep=jpk1597 DO jj = 2, jpjm11598 DO ji = 2, jpim11599 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )1600 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)1601 IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj ) - 1)1602 IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj ) - 1)1603 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji ,jj-1) - 1)1604 IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji ,jj+1) - 1)1605 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1606 IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN1607 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp1608 END IF1609 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN1610 misfdep(ji,jj) = ibtest1611 risfdep(ji,jj) = gdepw_1d(ibtest)1612 ENDIF1613 ENDDO1614 ENDDO1615 1616 IF( lk_mpp ) THEN1617 zbathy(:,:) = FLOAT( misfdep(:,:) )1618 CALL lbc_lnk( zbathy, 'T', 1. )1619 misfdep(:,:) = INT( zbathy(:,:) )1620 CALL lbc_lnk( risfdep, 'T', 1. )1621 CALL lbc_lnk( bathy, 'T', 1. )1622 zbathy(:,:) = FLOAT( mbathy(:,:) )1623 CALL lbc_lnk( zbathy, 'T', 1. )1624 mbathy(:,:) = INT( zbathy(:,:) )1625 ENDIF1626 !1627 !! fill hole in bathymetry1628 zmbathy (:,:)=mbathy (:,:)1629 DO jj = 2, jpjm11630 DO ji = 2, jpim11631 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )1632 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)1633 IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj ) + 1)1634 IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 01635 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 01636 IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 01637 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)1638 IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN1639 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;1640 END IF1641 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN1642 mbathy(ji,jj) = ibtest1643 bathy(ji,jj) = gdepw_1d(ibtest+1)1644 ENDIF1645 END DO1646 END DO1647 IF( lk_mpp ) THEN1648 zbathy(:,:) = FLOAT( misfdep(:,:) )1649 CALL lbc_lnk( zbathy, 'T', 1. )1650 misfdep(:,:) = INT( zbathy(:,:) )1651 CALL lbc_lnk( risfdep, 'T', 1. )1652 CALL lbc_lnk( bathy, 'T', 1. )1653 zbathy(:,:) = FLOAT( mbathy(:,:) )1654 CALL lbc_lnk( zbathy, 'T', 1. )1655 mbathy(:,:) = INT( zbathy(:,:) )1656 ENDIF1657 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1658 DO jj = 1, jpjm11659 DO ji = 1, jpim11660 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1661 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1662 END IF1663 END DO1664 END DO1665 IF( lk_mpp ) THEN1666 zbathy(:,:) = FLOAT( misfdep(:,:) )1667 CALL lbc_lnk( zbathy, 'T', 1. )1668 misfdep(:,:) = INT( zbathy(:,:) )1669 CALL lbc_lnk( risfdep, 'T', 1. )1670 CALL lbc_lnk( bathy, 'T', 1. )1671 zbathy(:,:) = FLOAT( mbathy(:,:) )1672 CALL lbc_lnk( zbathy, 'T', 1. )1673 mbathy(:,:) = INT( zbathy(:,:) )1674 ENDIF1675 ! if not compatible after all check (ie U point water column less than 2 cells), mask U1676 DO jj = 1, jpjm11677 DO ji = 1, jpim11678 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN1679 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;1680 END IF1681 END DO1682 END DO1683 IF( lk_mpp ) THEN1684 zbathy(:,:) = FLOAT( misfdep(:,:) )1685 CALL lbc_lnk( zbathy, 'T', 1. )1686 misfdep(:,:) = INT( zbathy(:,:) )1687 CALL lbc_lnk( risfdep, 'T', 1. )1688 CALL lbc_lnk( bathy, 'T', 1. )1689 zbathy(:,:) = FLOAT( mbathy(:,:) )1690 CALL lbc_lnk( zbathy, 'T', 1. )1691 mbathy(:,:) = INT( zbathy(:,:) )1692 ENDIF1693 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1694 DO jj = 1, jpjm11695 DO ji = 1, jpi1696 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1697 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;1698 END IF1699 END DO1700 END DO1701 IF( lk_mpp ) THEN1702 zbathy(:,:) = FLOAT( misfdep(:,:) )1703 CALL lbc_lnk( zbathy, 'T', 1. )1704 misfdep(:,:) = INT( zbathy(:,:) )1705 CALL lbc_lnk( risfdep, 'T', 1. )1706 CALL lbc_lnk( bathy, 'T', 1. )1707 zbathy(:,:) = FLOAT( mbathy(:,:) )1708 CALL lbc_lnk( zbathy, 'T', 1. )1709 mbathy(:,:) = INT( zbathy(:,:) )1710 ENDIF1711 ! if not compatible after all check (ie V point water column less than 2 cells), mask V1712 DO jj = 1, jpjm11713 DO ji = 1, jpi1714 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN1715 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;1716 END IF1717 END DO1718 END DO1719 IF( lk_mpp ) THEN1720 zbathy(:,:) = FLOAT( misfdep(:,:) )1721 CALL lbc_lnk( zbathy, 'T', 1. )1722 misfdep(:,:) = INT( zbathy(:,:) )1723 CALL lbc_lnk( risfdep, 'T', 1. )1724 CALL lbc_lnk( bathy, 'T', 1. )1725 zbathy(:,:) = FLOAT( mbathy(:,:) )1726 CALL lbc_lnk( zbathy, 'T', 1. )1727 mbathy(:,:) = INT( zbathy(:,:) )1728 ENDIF1729 ! if not compatible after all check, mask T1730 DO jj = 1, jpj1731 DO ji = 1, jpi1732 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN1733 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;1734 END IF1735 END DO1736 END DO1737 1738 WHERE (mbathy(:,:) == 1)1739 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp1740 END WHERE1741 END DO1742 ! end check compatibility ice shelf/bathy1743 ! remove very shallow ice shelf (less than ~ 10m if 75L)1744 WHERE (misfdep(:,:) <= 5)1745 misfdep = 1; risfdep = 0.0_wp;1746 END WHERE1747 1748 IF( icompt == 0 ) THEN1749 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'1750 ELSE1751 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'1752 ENDIF1753 1831 1754 1832 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1755 1833 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1756 1757 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1758 1759 END SUBROUTINE 1834 ! 1835 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1836 ! 1837 END SUBROUTINE zgr_isf 1838 1760 1839 1761 1840 SUBROUTINE zgr_sco … … 1803 1882 !! 1804 1883 !!---------------------------------------------------------------------- 1805 !1806 1884 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1807 1885 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers … … 1812 1890 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1813 1891 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1814 1892 !! 1815 1893 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1816 1894 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1817 1895 !!---------------------------------------------------------------------- 1818 1896 ! 1819 1897 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1820 1898 ! 1821 CALL wrk_alloc( jpi, jpj,zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )1899 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1822 1900 ! 1823 1901 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 1881 1959 DO jj = 1, jpj 1882 1960 DO ji = 1, jpi 1883 IF( bathy(ji,jj) == 0._wp ) THEN 1884 iip1 = MIN( ji+1, jpi ) 1885 ijp1 = MIN( jj+1, jpj ) 1886 iim1 = MAX( ji-1, 1 ) 1887 ijm1 = MAX( jj-1, 1 ) 1888 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1889 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1890 zenv(ji,jj) = rn_sbot_min 1891 ENDIF 1961 IF( bathy(ji,jj) == 0._wp ) THEN 1962 iip1 = MIN( ji+1, jpi ) 1963 ijp1 = MIN( jj+1, jpj ) 1964 iim1 = MAX( ji-1, 1 ) 1965 ijm1 = MAX( jj-1, 1 ) 1966 !!gm BUG fix see ticket #1617 1967 IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) & 1968 & + bathy(iim1,jj ) + bathy(iip1,jj ) & 1969 & + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) zenv(ji,jj) = rn_sbot_min 1970 !!gm 1971 !!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + & 1972 !!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1973 !!gm zenv(ji,jj) = rn_sbot_min 1974 !!gm ENDIF 1975 !!gm end 1892 1976 ENDIF 1893 1977 END DO … … 1983 2067 ENDIF 1984 2068 ! 1985 IF(lwp) THEN ! Control print1986 WRITE(numout,*)1987 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'1988 WRITE(numout,*)1989 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )1990 IF( nprint == 1 ) THEN1991 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )1992 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )1993 ENDIF1994 ENDIF1995 1996 2069 ! ! ============================== 1997 2070 ! ! hbatu, hbatv, hbatf fields … … 2111 2184 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 2112 2185 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2113 2114 fsdepw(:,:,:) = gdepw_0 (:,:,:)2115 fsde3w(:,:,:) = gdep3w_0(:,:,:)2116 2186 ! 2117 2187 IF( .NOT.ln_wd ) THEN 2118 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.02119 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.02120 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.02121 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.02122 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.02123 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.02124 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.02188 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2189 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2190 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2191 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2192 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2193 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2194 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2125 2195 END IF 2126 2196 2127 2197 #if defined key_agrif 2128 ! Ensure meaningful vertical scale factors in ghost lines/columns 2129 IF( .NOT. Agrif_Root() ) THEN 2130 ! 2131 IF((nbondi == -1).OR.(nbondi == 2)) THEN 2132 e3u_0(1,:,:) = e3u_0(2,:,:) 2133 ENDIF 2134 ! 2135 IF((nbondi == 1).OR.(nbondi == 2)) THEN 2136 e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) 2137 ENDIF 2138 ! 2139 IF((nbondj == -1).OR.(nbondj == 2)) THEN 2140 e3v_0(:,1,:) = e3v_0(:,2,:) 2141 ENDIF 2142 ! 2143 IF((nbondj == 1).OR.(nbondj == 2)) THEN 2144 e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) 2145 ENDIF 2146 ! 2147 ENDIF 2198 IF( .NOT. Agrif_Root() ) THEN ! Ensure meaningful vertical scale factors in ghost lines/columns 2199 IF( nbondi == -1 .OR. nbondi == 2 ) e3u_0( 1 , : ,:) = e3u_0( 2 , : ,:) 2200 IF( nbondi == 1 .OR. nbondi == 2 ) e3u_0(nlci-1, : ,:) = e3u_0(nlci-2, : ,:) 2201 IF( nbondj == -1 .OR. nbondj == 2 ) e3v_0( : , 1 ,:) = e3v_0( : , 2 ,:) 2202 IF( nbondj == 1 .OR. nbondj == 2 ) e3v_0( : ,nlcj-1,:) = e3v_0( : ,nlcj-2,:) 2203 ENDIF 2148 2204 #endif 2149 2205 2150 fsdept(:,:,:) = gdept_0 (:,:,:) 2151 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2152 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2153 fse3t (:,:,:) = e3t_0 (:,:,:) 2154 fse3u (:,:,:) = e3u_0 (:,:,:) 2155 fse3v (:,:,:) = e3v_0 (:,:,:) 2156 fse3f (:,:,:) = e3f_0 (:,:,:) 2157 fse3w (:,:,:) = e3w_0 (:,:,:) 2158 fse3uw(:,:,:) = e3uw_0 (:,:,:) 2159 fse3vw(:,:,:) = e3vw_0 (:,:,:) 2206 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 2207 !!gm and only that !!!!! 2208 !!gm THIS should be removed from here ! 2209 gdept_n(:,:,:) = gdept_0(:,:,:) 2210 gdepw_n(:,:,:) = gdepw_0(:,:,:) 2211 gde3w_n(:,:,:) = gde3w_0(:,:,:) 2212 e3t_n (:,:,:) = e3t_0 (:,:,:) 2213 e3u_n (:,:,:) = e3u_0 (:,:,:) 2214 e3v_n (:,:,:) = e3v_0 (:,:,:) 2215 e3f_n (:,:,:) = e3f_0 (:,:,:) 2216 e3w_n (:,:,:) = e3w_0 (:,:,:) 2217 e3uw_n (:,:,:) = e3uw_0 (:,:,:) 2218 e3vw_n (:,:,:) = e3vw_0 (:,:,:) 2219 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 2220 !! gm end 2160 2221 !! 2161 2222 ! HYBRID : … … 2163 2224 DO ji = 1, jpi 2164 2225 DO jk = 1, jpkm1 2165 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2166 END DO 2167 2226 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2227 END DO 2168 2228 IF( ln_wd ) THEN 2169 2229 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN … … 2184 2244 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2185 2245 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2186 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde p3w_0(:,:,:) )2187 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 2188 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 2189 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 2246 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 2247 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2248 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 2249 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 2190 2250 & ' w ', MINVAL( e3w_0 (:,:,:) ) 2191 2251 2192 2252 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2193 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde p3w_0(:,:,:) )2194 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 2195 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 2196 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 2253 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 2254 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2255 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 2256 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 2197 2257 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 2198 2258 ENDIF … … 2226 2286 END DO 2227 2287 ENDIF 2228 2229 !================================================================================2230 ! check the coordinate makes sense2231 !================================================================================2288 ! 2289 !================================================================================ 2290 ! check the coordinate makes sense 2291 !================================================================================ 2232 2292 DO ji = 1, jpi 2233 2293 DO jj = 1, jpj 2234 2294 ! 2235 2295 IF( hbatt(ji,jj) > 0._wp) THEN 2236 2296 DO jk = 1, mbathy(ji,jj) 2237 2297 ! check coordinate is monotonically increasing 2238 IF ( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN2298 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 2239 2299 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2240 2300 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2241 WRITE(numout,*) 'e3w', fse3w(ji,jj,:)2242 WRITE(numout,*) 'e3t', fse3t(ji,jj,:)2301 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2302 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 2243 2303 CALL ctl_stop( ctmp1 ) 2244 2304 ENDIF 2245 2305 ! and check it has never gone negative 2246 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN2306 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 2247 2307 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2248 2308 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2249 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2250 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2309 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2310 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2251 2311 CALL ctl_stop( ctmp1 ) 2252 2312 ENDIF 2253 2313 ! and check it never exceeds the total depth 2254 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN2314 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2255 2315 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2256 2316 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2257 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2317 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2258 2318 CALL ctl_stop( ctmp1 ) 2259 2319 ENDIF 2260 2320 END DO 2261 2321 ! 2262 2322 DO jk = 1, mbathy(ji,jj)-1 2263 2323 ! and check it never exceeds the total depth 2264 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN2324 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2265 2325 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2266 2326 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2267 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2327 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2268 2328 CALL ctl_stop( ctmp1 ) 2269 2329 ENDIF 2270 2330 END DO 2271 2272 2331 ENDIF 2273 2274 2332 END DO 2275 2333 END DO … … 2281 2339 END SUBROUTINE zgr_sco 2282 2340 2283 !!====================================================================== 2341 2284 2342 SUBROUTINE s_sh94() 2285 2286 2343 !!---------------------------------------------------------------------- 2287 2344 !! *** ROUTINE s_sh94 *** … … 2294 2351 !! Reference : Song and Haidvogel 1994. 2295 2352 !!---------------------------------------------------------------------- 2296 !2297 2353 INTEGER :: ji, jj, jk ! dummy loop argument 2298 2354 REAL(wp) :: zcoeft, zcoefw ! temporary scalars … … 2302 2358 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2303 2359 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2304 2305 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2306 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2360 !!---------------------------------------------------------------------- 2361 2362 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2363 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2307 2364 2308 2365 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp … … 2310 2367 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2311 2368 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2312 2369 ! 2313 2370 DO ji = 1, jpi 2314 2371 DO jj = 1, jpj 2315 2372 ! 2316 2373 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 2317 2374 DO jk = 1, jpk … … 2342 2399 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2343 2400 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2344 gdept_0 2345 gdepw_0 2346 gde p3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )2401 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2402 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2403 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 2347 2404 END DO 2348 2405 ! … … 2403 2460 END DO 2404 2461 END DO 2405 2406 CALL wrk_dealloc( jpi, jpj, jpk,z_gsigw3, z_gsigt3, z_gsi3w3 )2407 CALL wrk_dealloc( jpi, jpj, jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )2408 2462 ! 2463 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2464 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2465 ! 2409 2466 END SUBROUTINE s_sh94 2410 2467 2468 2411 2469 SUBROUTINE s_sf12 2412 2413 2470 !!---------------------------------------------------------------------- 2414 2471 !! *** ROUTINE s_sf12 *** … … 2426 2483 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 2427 2484 !!---------------------------------------------------------------------- 2428 !2429 2485 INTEGER :: ji, jj, jk ! dummy loop argument 2430 2486 REAL(wp) :: zsmth ! smoothing around critical depth … … 2435 2491 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2436 2492 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2437 2493 !!---------------------------------------------------------------------- 2438 2494 ! 2439 2495 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) … … 2499 2555 2500 2556 DO jk = 1, jpk 2501 gdept_0 2502 gdepw_0 2503 gde p3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)2557 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2558 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2559 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 2504 2560 END DO 2505 2561 … … 2573 2629 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 2574 2630 ! 2575 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)2631 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 2576 2632 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 2577 2633 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) … … 2586 2642 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 2587 2643 ! 2588 ! ! ============= 2589 2590 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2591 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2592 2644 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2645 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2646 ! 2593 2647 END SUBROUTINE s_sf12 2594 2648 2649 2595 2650 SUBROUTINE s_tanh() 2596 2597 2651 !!---------------------------------------------------------------------- 2598 2652 !! *** ROUTINE s_tanh*** … … 2604 2658 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 2605 2659 !!---------------------------------------------------------------------- 2606 2607 INTEGER :: ji, jj, jk ! dummy loop argument 2660 INTEGER :: ji, jj, jk ! dummy loop argument 2608 2661 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2609 2610 2662 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 2611 2663 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 2612 2613 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2614 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2664 !!---------------------------------------------------------------------- 2665 2666 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2667 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2615 2668 2616 2669 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp … … 2642 2695 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2643 2696 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2644 gdept_0 2645 gdepw_0 2646 gde p3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )2697 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2698 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2699 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2647 2700 END DO 2648 2701 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 2661 2714 END DO 2662 2715 END DO 2663 2664 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w)2665 CALL wrk_dealloc( jpk, z_esigt, z_esigw)2666 2716 ! 2717 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2718 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2719 ! 2667 2720 END SUBROUTINE s_tanh 2721 2668 2722 2669 2723 FUNCTION fssig( pk ) RESULT( pf ) … … 2737 2791 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 2738 2792 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 2739 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth2740 REAL(wp), INTENT(in ) :: pzs ! surface box depth2741 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter2742 REAL(wp) :: za1,za2,za3 ! local variables2743 REAL(wp) :: zn1,zn2 ! local variables2744 REAL(wp) :: za,zb,zx ! local variables2745 integer :: jk2746 !!----------------------------------------------------------------------2747 ! 2748 2749 zn1 = 1. /(jpk-1.)2750 zn2 = 1. - zn12751 2793 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 2794 REAL(wp), INTENT(in ) :: pzs ! surface box depth 2795 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 2796 ! 2797 INTEGER :: jk ! dummy loop index 2798 REAL(wp) :: za1,za2,za3 ! local scalar 2799 REAL(wp) :: zn1,zn2 ! - - 2800 REAL(wp) :: za,zb,zx ! - - 2801 !!---------------------------------------------------------------------- 2802 ! 2803 zn1 = 1._wp / REAL( jpkm1, wp ) 2804 zn2 = 1._wp - zn1 2805 ! 2752 2806 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 2753 2807 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 2754 2808 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2755 2809 ! 2756 2810 za = pzb - za3*(pzs-za1)-za2 2757 2811 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 2758 2812 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 2759 2813 zx = 1.0_wp-za/2.0_wp-zb 2760 2814 ! 2761 2815 DO jk = 1, jpk 2762 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + &2763 &zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &2764 &(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )2816 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 2817 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 2818 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 2765 2819 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 2766 ENDDO 2767 2820 END DO 2768 2821 ! 2769 2822 END FUNCTION fgamma -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r5870 r6141 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_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r6138 r6141 35 35 USE dtauvd ! data: U & V current (dta_uvd routine) 36 36 USE domvvl ! varying vertical mesh 37 USE iscplrst ! ice sheet coupling 37 38 ! 38 39 USE in_out_manager ! I/O manager … … 49 50 50 51 !! * Substitutions 51 # include "domzgr_substitute.h90"52 52 # include "vectopt_loop_substitute.h90" 53 53 !!---------------------------------------------------------------------- … … 85 85 IF( ln_rstart ) THEN ! Restart from a file 86 86 ! ! ------------------- 87 CALL rst_read ! Read the restart file 88 CALL day_init ! model calendar (using both namelist and restart infos) 87 CALL rst_read ! Read the restart file 88 IF (ln_iscpl) CALL iscpl_stp ! extraloate restart to wet and dry 89 CALL day_init ! model calendar (using both namelist and restart infos) 89 90 ELSE 90 91 ! ! Start from rest … … 120 121 ENDIF 121 122 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 123 ! 124 !!gm This is to be changed !!!! 125 ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here 126 IF( .NOT.ln_linssh ) THEN 125 127 DO jk = 1, jpk 126 fse3t_b(:,:,jk) = fse3t_n(:,:,jk)128 e3t_b(:,:,jk) = e3t_n(:,:,jk) 127 129 END DO 128 130 ENDIF 131 !!gm 129 132 ! 130 133 ENDIF 131 !132 134 ! 133 135 ! Initialize "now" and "before" barotropic velocities: 134 ! Do it whatever the free surface method, these arrays 135 ! being eventually used 136 ! 136 ! Do it whatever the free surface method, these arrays being eventually used 137 137 ! 138 138 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 139 139 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 140 140 ! 141 !!gm the use of umsak & vmask is not necessary belox as un, vn, ub, vb are always masked 141 142 DO jk = 1, jpkm1 142 143 DO jj = 1, jpj 143 144 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)145 un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 146 vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 146 147 ! 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)148 ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 149 vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 149 150 END DO 150 151 END DO 151 152 END DO 152 153 ! 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 ! 154 un_b(:,:) = un_b(:,:) * r1_hu_n(:,:) 155 vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:) 156 ! 157 ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 158 vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 159 159 ! 160 160 IF( nn_timing == 1 ) CALL timing_stop('istate_init') … … 184 184 ! 185 185 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)186 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((gdept_n(:,:,jk)-80.)/30.) ) & 187 & + 10. * ( 5000. - gdept_n(:,:,jk) ) /5000.) ) * tmask(:,:,jk) 188 188 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 189 189 END DO … … 238 238 ! 239 239 DO jk = 1, jpk 240 tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)240 tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - gdept_n(:,:,jk) / 1000 ) ) * tmask(:,:,jk) 241 241 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 242 242 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 243 ! 248 244 ! set salinity field to a constant value … … 314 310 tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb 315 311 ! 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 312 ! set salinity field to a constant value 321 313 ! -------------------------------------- … … 363 355 DO jj = 1, jpj 364 356 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) / 2357 tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 ) ) & 358 & * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2 & 359 & + ( 15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) ) & 360 & - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.) & 361 & + 7. * (1500. - gdept_n(ji,jj,jk)) / 1500. ) & 362 & * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 371 363 tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 372 364 tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) 373 365 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) / 2366 tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 ) ) & 367 & * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2 & 368 & + ( 35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000. & 369 & - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60. ) / 650. ) & 370 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 35. ) / 100. ) & 371 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.) ) & 372 & * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 381 373 tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 382 374 tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) … … 451 443 zalfg = 0.5 * grav * rau0 452 444 453 zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value445 zprn(:,:,1) = zalfg * e3w_n(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value 454 446 455 447 DO jk = 2, jpkm1 ! Vertical integration from the surface 456 448 zprn(:,:,jk) = zprn(:,:,jk-1) & 457 & + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )449 & + zalfg * e3w_n(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) ) 458 450 END DO 459 451
Note: See TracChangeset
for help on using the changeset viewer.