Changeset 9939 for NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM
- Timestamp:
- 2018-07-13T09:28:50+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3
- Files:
-
- 8 edited
- 2 copied
-
. (copied) (copied from NEMO/trunk)
-
src/OCE/DOM/daymod.F90 (modified) (6 diffs)
-
src/OCE/DOM/dom_oce.F90 (modified) (2 diffs)
-
src/OCE/DOM/domain.F90 (modified) (5 diffs)
-
src/OCE/DOM/domvvl.F90 (modified) (33 diffs)
-
src/OCE/DOM/iscplini.F90 (modified) (2 diffs)
-
src/OCE/DOM/iscplrst.F90 (modified) (1 diff)
-
src/OCE/DOM/istate.F90 (modified) (1 diff)
-
src/OCE/DOM/phycst.F90 (modified) (3 diffs)
-
src/OCE/DOM/restart.F90 (copied) (copied from NEMO/trunk/src/OCE/IOM/restart.F90) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/daymod.F90
r9598 r9939 20 20 !! ------------------------------- 21 21 !! sbcmod assume that the time step is dividing the number of second of 22 !! in a day, i.e. ===> MOD( rday, r dt ) == 022 !! in a day, i.e. ===> MOD( rday, rn_Dt ) == 0 23 23 !! except when user defined forcing is used (see sbcmod.F90) 24 24 !!---------------------------------------------------------------------- … … 72 72 ! 73 73 ! max number of seconds between each restart 74 IF( REAL( nitend - nit000 + 1 ) * r dt > REAL( HUGE( nsec1jan000 ) ) ) THEN74 IF( REAL( nitend - nit000 + 1 ) * rn_Dt > REAL( HUGE( nsec1jan000 ) ) ) THEN 75 75 CALL ctl_stop( 'The number of seconds between each restart exceeds the integer 4 max value: 2^31-1. ', & 76 76 & 'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' ) 77 77 ENDIF 78 nsecd = NINT( rday )79 nsecd05 = NINT( 0.5 * rday )80 ndt = NINT( r dt)81 ndt05 = NINT( 0.5 * r dt)78 nsecd = NINT( rday ) 79 nsecd05 = NINT( 0.5 * rday ) 80 ndt = NINT( rn_Dt ) 81 ndt05 = NINT( 0.5 * rn_Dt ) 82 82 83 83 IF( .NOT. l_offline ) CALL day_rst( nit000, 'READ' ) … … 239 239 nsec_week = nsec_week + ndt 240 240 nsec_day = nsec_day + ndt 241 adatrj = adatrj + r dt / rday242 fjulday = fjulday + r dt / rday241 adatrj = adatrj + rn_Dt / rday 242 fjulday = fjulday + rn_Dt / rday 243 243 IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < zprec ) fjulday = REAL(NINT(fjulday),wp) ! avoid truncation error 244 244 IF( ABS(adatrj - REAL(NINT(adatrj ),wp)) < zprec ) adatrj = REAL(NINT(adatrj ),wp) ! avoid truncation error … … 309 309 !! In both those options, the exact duration of the experiment 310 310 !! since the beginning (cumulated duration of all previous restart runs) 311 !! is not stored in the restart and is assumed to be (nit000-1)*r dt.311 !! is not stored in the restart and is assumed to be (nit000-1)*rn_Dt. 312 312 !! This is valid is the time step has remained constant. 313 313 !! … … 378 378 nminute = ( nn_time0 - nhour * 100 ) 379 379 IF( nhour*3600+nminute*60-ndt05 .lt. 0 ) ndastp=ndastp-1 ! Start hour is specified in the namelist (default 0) 380 adatrj = ( REAL( nit000-1, wp ) * r dt ) / rday380 adatrj = ( REAL( nit000-1, wp ) * rn_Dt ) / rday 381 381 ! note this is wrong if time step has changed during run 382 382 ENDIF … … 387 387 nminute = ( nn_time0 - nhour * 100 ) 388 388 IF( nhour*3600+nminute*60-ndt05 .lt. 0 ) ndastp=ndastp-1 ! Start hour is specified in the namelist (default 0) 389 adatrj = ( REAL( nit000-1, wp ) * r dt ) / rday389 adatrj = ( REAL( nit000-1, wp ) * rn_Dt ) / rday 390 390 ENDIF 391 391 IF( ABS(adatrj - REAL(NINT(adatrj),wp)) < 0.1 / rday ) adatrj = REAL(NINT(adatrj),wp) ! avoid truncation error -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90
r9667 r9939 33 33 LOGICAL , PUBLIC :: ln_meshmask !: =T create a mesh-mask file (mesh_mask.nc) 34 34 REAL(wp), PUBLIC :: rn_isfhmin !: threshold to discriminate grounded ice to floating ice 35 REAL(wp), PUBLIC :: rn_ rdt!: time step for the dynamics and tracer35 REAL(wp), PUBLIC :: rn_dt !: time step for the dynamics and tracer 36 36 REAL(wp), PUBLIC :: rn_atfp !: asselin time filter parameter 37 INTEGER , PUBLIC :: nn_euler!: =0 start with forward time step or not (=1)37 LOGICAL , PUBLIC :: ln_1st_euler !: =0 start with forward time step or not (=1) 38 38 LOGICAL , PUBLIC :: ln_iscpl !: coupling with ice sheet 39 39 LOGICAL , PUBLIC :: ln_crs !: Apply grid coarsening to dynamical model output or online passive tracers … … 50 50 LOGICAL, PUBLIC :: ln_bt_auto !: Set number of barotropic iterations automatically 51 51 INTEGER, PUBLIC :: nn_bt_flt !: Filter choice 52 INTEGER, PUBLIC :: nn_ baro !: Number of barotropic iterations during one baroclinic step (rdt)52 INTEGER, PUBLIC :: nn_e !: Number of external mode sub-step used at each ocean time-step 53 53 REAL(wp), PUBLIC :: rn_bt_cmax !: Maximum allowed courant number (used if ln_bt_auto=T) 54 54 REAL(wp), PUBLIC :: rn_bt_alpha !: Time stepping diffusion parameter 55 55 56 57 ! !! old non-DOCTOR names still used in the model58 REAL(wp), PUBLIC :: atfp !: asselin time filter parameter59 REAL(wp), PUBLIC :: rdt !: time step for the dynamics and tracer60 61 56 ! !!! associated variables 62 INTEGER , PUBLIC :: neuler !: restart euler forward option (0=Euler) 63 REAL(wp), PUBLIC :: r2dt !: = 2*rdt except at nit000 (=rdt) if neuler=0 57 LOGICAL , PUBLIC :: l_1st_euler !: Euler 1st time-step flag (=T if ln_restart=F or ln_1st_euler=T) 58 REAL(wp), PUBLIC :: rDt, r1_Dt !: MLF: = 2*rn_Dt and 1/(2*rn_Dt) except if l_1st_euler=T where half the value is used 59 ! ! RK3: = rn_Dt 64 60 65 61 !!---------------------------------------------------------------------- -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90
r9598 r9939 288 288 INTEGER :: ios ! Local integer 289 289 ! 290 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, &291 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl , &292 & nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate , &293 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, nn_euler, &290 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, & 291 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl , & 292 & nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate , & 293 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, ln_1st_euler, & 294 294 & ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 295 NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_ rdt, rn_atfp, ln_crs, ln_meshmask295 NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_Dt, rn_atfp, ln_crs, ln_meshmask 296 296 #if defined key_netcdf4 297 297 NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip … … 323 323 WRITE(numout,*) ' restart output directory cn_ocerst_outdir= ', TRIM( cn_ocerst_outdir ) 324 324 WRITE(numout,*) ' restart logical ln_rstart = ', ln_rstart 325 WRITE(numout,*) ' start with forward time step nn_euler = ', nn_euler325 WRITE(numout,*) ' start with forward time step ln_1st_euler = ', ln_1st_euler 326 326 WRITE(numout,*) ' control of time step nn_rstctl = ', nn_rstctl 327 327 WRITE(numout,*) ' number of the first time step nn_it000 = ', nn_it000 … … 361 361 nstocklist = nn_stocklist 362 362 nwrite = nn_write 363 neuler = nn_euler 364 IF( neuler == 1 .AND. .NOT. ln_rstart ) THEN 363 IF( ln_rstart ) THEN ! restart : set 1st time-step scheme (LF or forward) 364 l_1st_euler = ln_1st_euler 365 ELSE ! start from rest : always an Euler scheme for the 1st time-step 365 366 IF(lwp) WRITE(numout,*) 366 367 IF(lwp) WRITE(numout,*)' ==>>> Start from rest (ln_rstart=F)' 367 IF(lwp) WRITE(numout,*)' an Euler initial time step is used : nn_euler is forced to 0'368 neuler = 0368 IF(lwp) WRITE(numout,*)' an Euler initial time step is used ' 369 l_1st_euler = .TRUE. 369 370 ENDIF 370 371 ! ! control of output frequency … … 374 375 nstock = nitend 375 376 ENDIF 376 IF 377 IF( nwrite == 0 ) THEN 377 378 WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend 378 379 CALL ctl_warn( ctmp1 ) … … 413 414 WRITE(numout,*) ' create mesh/mask file ln_meshmask = ', ln_meshmask 414 415 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' [m]' 415 WRITE(numout,*) ' ocean time step rn_ rdt = ', rn_rdt416 WRITE(numout,*) ' ocean time step rn_dt = ', rn_dt 416 417 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 417 418 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs 418 419 ENDIF 419 420 ! 420 ! ! conversion DOCTOR names into model names (this should disappear soon)421 atfp = rn_atfp422 rdt = rn_rdt423 424 421 IF( TRIM(Agrif_CFixed()) == '0' ) THEN 425 422 lrxios = ln_xios_read.AND.ln_rstart -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90
r9598 r9939 54 54 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 55 55 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td 57 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf 58 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors59 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors60 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t 61 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 57 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 58 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n ! baroclinic scale factors 59 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_a ! baroclinic scale factors 60 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 61 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 62 62 63 63 !! * Substitutions … … 76 76 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 77 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 78 ALLOCATE( t ilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &79 & dt ilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &78 ALLOCATE( te3t_b(jpi,jpj,jpk) , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) , & 79 & dte3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , & 80 80 & STAT = dom_vvl_alloc ) 81 81 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) … … 103 103 !! - interpolate scale factors 104 104 !! 105 !! ** Action : - e3t_(n/b) and t ilde_e3t_(n/b)105 !! ** Action : - e3t_(n/b) and te3t_(n/b) 106 106 !! - Regrid: e3(u/v)_n 107 107 !! e3(u/v)_b … … 117 117 INTEGER :: ji, jj, jk 118 118 INTEGER :: ii0, ii1, ij0, ij1 119 REAL(wp):: zcoef 119 REAL(wp):: zcoef, z1_Dt 120 120 !!---------------------------------------------------------------------- 121 121 ! … … 129 129 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 130 130 ! 131 ! ! Read or initialize e3t_(b/n), t ilde_e3t_(b/n) and hdiv_lf131 ! ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf 132 132 CALL dom_vvl_rst( nit000, 'READ' ) 133 133 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all … … 208 208 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 209 209 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 210 frq_rst_hdv(:,:) = 1._wp / r dt210 frq_rst_hdv(:,:) = 1._wp / rn_Dt 211 211 ENDIF 212 212 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 213 z1_Dt = 1._wp / rn_Dt 213 214 DO jj = 1, jpj 214 215 DO ji = 1, jpi … … 216 217 IF( ABS(gphit(ji,jj)) >= 6.) THEN 217 218 ! values outside the equatorial band and transition zone (ztilde) 218 frq_rst_e3t(ji,jj) = 2. 0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp )219 frq_rst_hdv(ji,jj) = 2. 0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )219 frq_rst_e3t(ji,jj) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400._wp ) 220 frq_rst_hdv(ji,jj) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400._wp ) 220 221 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 221 222 ! values inside the equatorial band (ztilde as zstar) 222 frq_rst_e3t(ji,jj) = 0. 0_wp223 frq_rst_hdv(ji,jj) = 1.0_wp / rdt223 frq_rst_e3t(ji,jj) = 0._wp 224 frq_rst_hdv(ji,jj) = z1_Dt 224 225 ELSE ! transition band (2.5 to 6 degrees N/S) 225 226 ! ! (linearly transition from z-tilde to z-star) 226 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 227 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 228 & * 180._wp / 3.5_wp ) ) 229 frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 230 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 231 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 232 & * 180._wp / 3.5_wp ) ) 227 frq_rst_e3t(ji,jj) = 0._wp + ( frq_rst_e3t(ji,jj) - 0._wp ) * 0.5_wp & 228 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp ) ) 229 frq_rst_hdv(ji,jj) = z1_Dt + ( frq_rst_hdv(ji,jj) - z1_Dt ) * 0.5_wp & 230 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp ) ) 233 231 ENDIF 234 232 END DO … … 237 235 ii0 = 103 ; ii1 = 111 238 236 ij0 = 128 ; ij1 = 135 ; 239 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0. 0_wp240 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt237 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp 238 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = z1_Dt 241 239 ENDIF 242 240 ENDIF … … 280 278 !! 281 279 !! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case 282 !! - t ilde_e3t_a: after increment of vertical scale factor280 !! - te3t_a: after increment of vertical scale factor 283 281 !! in z_tilde case 284 282 !! - e3(t/u/v)_a … … 345 343 IF( kt > nit000 ) THEN 346 344 DO jk = 1, jpkm1 347 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - r dt * frq_rst_hdv(:,:) &345 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:) & 348 346 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 349 347 END DO … … 353 351 ! II - after z_tilde increments of vertical scale factors 354 352 ! ======================================================= 355 t ilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms353 te3t_a(:,:,:) = 0._wp ! te3t_a used to store tendency terms 356 354 357 355 ! 1 - High frequency divergence term … … 359 357 IF( ln_vvl_ztilde ) THEN ! z_tilde case 360 358 DO jk = 1, jpkm1 361 t ilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )359 te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 362 360 END DO 363 361 ELSE ! layer case 364 362 DO jk = 1, jpkm1 365 t ilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)363 te3t_a(:,:,jk) = te3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 366 364 END DO 367 365 ENDIF … … 371 369 IF( ln_vvl_ztilde ) THEN 372 370 DO jk = 1, jpk 373 t ilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)371 te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk) 374 372 END DO 375 373 ENDIF … … 383 381 DO ji = 1, fs_jpim1 ! vector opt. 384 382 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 385 & * ( t ilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )383 & * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj ,jk) ) 386 384 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 387 & * ( t ilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )385 & * ( te3t_b(ji,jj,jk) - te3t_b(ji ,jj+1,jk) ) 388 386 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 389 387 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) … … 400 398 DO jj = 2, jpjm1 401 399 DO ji = fs_2, fs_jpim1 ! vector opt. 402 t ilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) &400 te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 403 401 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 404 402 & ) * r1_e1e2t(ji,jj) … … 414 412 ! Leapfrog time stepping 415 413 ! ~~~~~~~~~~~~~~~~~~~~~~ 416 IF( neuler == 0 .AND. kt == nit000 ) THEN 417 z2dt = rdt 418 ELSE 419 z2dt = 2.0_wp * rdt 420 ENDIF 421 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 422 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 414 CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp ) 415 te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:) 423 416 424 417 ! Maximum deformation control … … 426 419 ze3t(:,:,jpk) = 0._wp 427 420 DO jk = 1, jpkm1 428 ze3t(:,:,jk) = t ilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)421 ze3t(:,:,jk) = te3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 429 422 END DO 430 423 z_tmax = MAXVAL( ze3t(:,:,:) ) … … 446 439 ENDIF 447 440 IF (lwp) THEN 448 WRITE(numout, *) 'MAX( t ilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax441 WRITE(numout, *) 'MAX( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 449 442 WRITE(numout, *) 'at i, j, k=', ijk_max 450 WRITE(numout, *) 'MIN( t ilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin443 WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 451 444 WRITE(numout, *) 'at i, j, k=', ijk_min 452 CALL ctl_warn('MAX( ABS( t ilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')445 CALL ctl_warn('MAX( ABS( te3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 453 446 ENDIF 454 447 ENDIF 455 448 ! - ML - end test 456 449 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 457 t ilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) )458 t ilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )450 te3t_a(:,:,:) = MIN( te3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) ) 451 te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 459 452 460 453 ! … … 462 455 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 463 456 DO jk = 1, jpkm1 464 dt ilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)457 dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_b(:,:,jk) 465 458 END DO 466 459 ! III - Barotropic repartition of the sea surface height over the baroclinic profile … … 470 463 ! i.e. locally and not spread over the water column. 471 464 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 472 zht(:,:) = 0. 465 zht(:,:) = 0._wp 473 466 DO jk = 1, jpkm1 474 zht(:,:) = zht(:,:) + t ilde_e3t_a(:,:,jk) * tmask(:,:,jk)467 zht(:,:) = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk) 475 468 END DO 476 469 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 477 470 DO jk = 1, jpkm1 478 dt ilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)471 dte3t_a(:,:,jk) = dte3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 479 472 END DO 480 473 … … 484 477 ! ! ---baroclinic part--------- ! 485 478 DO jk = 1, jpkm1 486 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dt ilde_e3t_a(:,:,jk) * tmask(:,:,jk)479 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dte3t_a(:,:,jk) * tmask(:,:,jk) 487 480 END DO 488 481 ENDIF … … 494 487 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 495 488 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 496 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(t ilde_e3t_a))) =', z_tmax489 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(te3t_a))) =', z_tmax 497 490 END IF 498 491 ! … … 573 566 !! - recompute depths and water height fields 574 567 !! 575 !! ** Action : - e3t_(b/n), t ilde_e3t_(b/n) and e3(u/v)_n ready for next time step568 !! ** Action : - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step 576 569 !! - Recompute: 577 570 !! e3(u/v)_b … … 587 580 INTEGER, INTENT( in ) :: kt ! time step 588 581 ! 589 INTEGER :: ji, jj, jk ! dummy loop indices590 REAL(wp) :: zcoef 582 INTEGER :: ji, jj, jk ! dummy loop indices 583 REAL(wp) :: zcoef, ze3f ! local scalar 591 584 !!---------------------------------------------------------------------- 592 585 ! … … 605 598 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 606 599 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 607 IF( neuler == 0 .AND. kt == nit000) THEN608 t ilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)600 IF( l_1st_euler ) THEN 601 te3t_n(:,:,:) = te3t_a(:,:,:) 609 602 ELSE 610 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 611 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 603 DO jk = 1, jpk 604 DO jj = 1, jpj 605 DO ji = 1, jpi 606 ze3f = te3t_n(ji,jj,jk) & 607 & + rn_atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) ) 608 te3t_b(ji,jj,jk) = ze3f 609 te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk) 610 END DO 611 END DO 612 END DO 612 613 ENDIF 613 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)614 614 ENDIF 615 615 gdept_b(:,:,:) = gdept_n(:,:,:) … … 806 806 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios ) 807 807 ! 808 id1 = iom_varid( numror, 'e3t_b' , ldstop = .FALSE. )809 id2 = iom_varid( numror, 'e3t_n' , ldstop = .FALSE. )808 id1 = iom_varid( numror, 'e3t_b' , ldstop = .FALSE. ) 809 id2 = iom_varid( numror, 'e3t_n' , ldstop = .FALSE. ) 810 810 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 811 811 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 812 id5 = iom_varid( numror, 'hdiv_lf' , ldstop = .FALSE. )812 id5 = iom_varid( numror, 'hdiv_lf' , ldstop = .FALSE. ) 813 813 ! ! --------- ! 814 814 ! ! all cases ! … … 823 823 e3t_b(:,:,:) = e3t_0(:,:,:) 824 824 END WHERE 825 IF( neuler == 0) THEN825 IF( l_1st_euler ) THEN 826 826 e3t_b(:,:,:) = e3t_n(:,:,:) 827 827 ENDIF … … 829 829 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 830 830 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 831 IF(lwp) write(numout,*) ' neuler is forced to 0'831 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 832 832 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 833 833 e3t_n(:,:,:) = e3t_b(:,:,:) 834 neuler = 0834 l_1st_euler = .TRUE. 835 835 ELSE IF( id2 > 0 ) THEN 836 836 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 837 837 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 838 IF(lwp) write(numout,*) ' neuler is forced to 0'838 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 839 839 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 840 840 e3t_b(:,:,:) = e3t_n(:,:,:) 841 neuler = 0841 l_1st_euler = .TRUE. 842 842 ELSE 843 843 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 844 844 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 845 IF(lwp) write(numout,*) ' neuler is forced to 0'845 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 846 846 DO jk = 1, jpk 847 847 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & … … 850 850 END DO 851 851 e3t_b(:,:,:) = e3t_n(:,:,:) 852 neuler = 0852 l_1st_euler = .TRUE. 853 853 ENDIF 854 854 ! ! ----------- ! … … 862 862 ! ! ----------------------- ! 863 863 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist 864 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', t ilde_e3t_b(:,:,:), ldxios = lrxios )865 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', t ilde_e3t_n(:,:,:), ldxios = lrxios )864 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lrxios ) 865 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios ) 866 866 ELSE ! one at least array is missing 867 t ilde_e3t_b(:,:,:) = 0.0_wp868 t ilde_e3t_n(:,:,:) = 0.0_wp867 te3t_b(:,:,:) = 0.0_wp 868 te3t_n(:,:,:) = 0.0_wp 869 869 ENDIF 870 870 ! ! ------------ ! … … 942 942 943 943 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 944 t ilde_e3t_b(:,:,:) = 0._wp945 t ilde_e3t_n(:,:,:) = 0._wp944 te3t_b(:,:,:) = 0._wp 945 te3t_n(:,:,:) = 0._wp 946 946 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp 947 947 END IF … … 960 960 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! 961 961 ! ! ----------------------- ! 962 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', t ilde_e3t_b(:,:,:), ldxios = lwxios)963 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', t ilde_e3t_n(:,:,:), ldxios = lwxios)962 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios) 963 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios) 964 964 END IF 965 965 ! ! -------------! … … 1016 1016 WRITE(numout,*) ' rn_rst_e3t = 0.e0' 1017 1017 WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 1018 WRITE(numout,*) ' rn_lf_cutoff = 1 .0/rdt'1018 WRITE(numout,*) ' rn_lf_cutoff = 1/rn_Dt' 1019 1019 ELSE 1020 1020 WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplini.F90
r9598 r9939 71 71 ! 72 72 nstp_iscpl=MIN( nn_fiscpl, nitend-nit000+1 ) ! the coupling period have to be less or egal than the total number of time step 73 rdt_iscpl = nstp_iscpl * rn_ rdt73 rdt_iscpl = nstp_iscpl * rn_Dt 74 74 ! 75 75 IF (lwp) THEN … … 79 79 WRITE(numout,*) ' conservation flag (ln_hsb ) = ', ln_hsb 80 80 WRITE(numout,*) ' nb of stp for cons (rn_fiscpl) = ', nstp_iscpl 81 IF (nstp_iscpl .NE.nn_fiscpl) WRITE(numout,*) 'W A R N I N G: nb of stp for cons has been modified &81 IF (nstp_iscpl /= nn_fiscpl) WRITE(numout,*) 'W A R N I N G: nb of stp for cons has been modified & 82 82 & (larger than run length)' 83 83 WRITE(numout,*) ' coupling time step = ', rdt_iscpl -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90
r9598 r9939 89 89 END IF 90 90 ! 91 neuler = 0! next step is an euler time step91 l_1st_euler = .TRUE. ! next step is an euler time step 92 92 ! 93 93 ! ! set _b and _n variables equal -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/istate.F90
r9598 r9939 92 92 ! ! --------------- 93 93 numror = 0 ! define numror = 0 -> no restart file to read 94 neuler = 0 ! Set time-step indicator at nit000 (euler forward)94 l_1st_euler = .TRUE. ! Set a Euler forward 1sr time-step 95 95 CALL day_init ! model calendar (using both namelist and restart infos) 96 96 ! ! Initialization of ocean to zero -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/phycst.F90
r9656 r9939 34 34 REAL(wp), PUBLIC :: rhhmm = 60._wp !: number of minutes in one hour 35 35 REAL(wp), PUBLIC :: rmmss = 60._wp !: number of seconds in one minute 36 REAL(wp), PUBLIC :: omega !: earth rotation parameter [s-1]37 REAL(wp), PUBLIC :: ra = 6371229._wp !: earth radius [m]38 REAL(wp), PUBLIC :: grav = 9.80665_wp !: gravity [m/s2]36 REAL(wp), PUBLIC :: omega !: earth rotation parameter [s-1] 37 REAL(wp), PUBLIC :: ra = 6371229._wp !: earth radius [m] 38 REAL(wp), PUBLIC :: grav = 9.80665_wp !: gravity [m/s2] 39 39 40 REAL(wp), PUBLIC :: rtt = 273.16_wp !: triple point of temperature [Kelvin] 41 REAL(wp), PUBLIC :: rt0 = 273.15_wp !: freezing point of fresh water [Kelvin] 42 REAL(wp), PUBLIC :: rt0_snow = 273.15_wp !: melting point of snow [Kelvin] 43 #if defined key_si3 44 REAL(wp), PUBLIC :: rt0_ice = 273.15_wp !: melting point of ice [Kelvin] 45 #else 46 REAL(wp), PUBLIC :: rt0_ice = 273.05_wp !: melting point of ice [Kelvin] 47 #endif 48 REAL(wp), PUBLIC :: rau0 !: volumic mass of reference [kg/m3] 49 REAL(wp), PUBLIC :: r1_rau0 !: = 1. / rau0 [m3/kg] 50 REAL(wp), PUBLIC :: rcp !: ocean specific heat [J/Kelvin] 51 REAL(wp), PUBLIC :: r1_rcp !: = 1. / rcp [Kelvin/J] 52 REAL(wp), PUBLIC :: rau0_rcp !: = rau0 * rcp 53 REAL(wp), PUBLIC :: r1_rau0_rcp !: = 1. / ( rau0 * rcp ) 40 REAL(wp), PUBLIC :: rt0 = 273.15_wp !: freezing point of fresh water [Kelvin] 41 REAL(wp), PUBLIC :: rho0 !: volumic mass of reference [kg/m3] 42 REAL(wp), PUBLIC :: r1_rho0 !: = 1. / rho0 [m3/kg] 43 REAL(wp), PUBLIC :: rcp !: ocean specific heat [J/Kelvin] 44 REAL(wp), PUBLIC :: r1_rcp !: = 1. / rcp [Kelvin/J] 45 REAL(wp), PUBLIC :: rho0_rcp !: = rho0 * rcp 46 REAL(wp), PUBLIC :: r1_rho0_rcp !: = 1. / ( rho0 * rcp ) 54 47 55 REAL(wp), PUBLIC :: rhosn = 330._wp !: volumic mass of snow [kg/m3] 56 REAL(wp), PUBLIC :: rhofw = 1000._wp !: volumic mass of freshwater in melt ponds [kg/m3] 48 REAL(wp), PUBLIC :: rhoi = 917._wp !: sea ice density [kg/m3] 49 REAL(wp), PUBLIC :: rhos = 330._wp !: snow density [kg/m3] 50 REAL(wp), PUBLIC :: rhow = 1000._wp !: water density (in melt ponds) [kg/m3] 51 REAL(wp), PUBLIC :: rcnd_i = 2.034396_wp !: thermal conductivity of fresh ice [W/m/K] 52 REAL(wp), PUBLIC :: rcpi = 2067.0_wp !: specific heat of fresh ice [J/kg/K] 53 REAL(wp), PUBLIC :: rLsub = 2.834e+6_wp !: pure ice latent heat of sublimation [J/kg] 54 REAL(wp), PUBLIC :: rLfus = 0.334e+6_wp !: latent heat of fusion of fresh ice [J/kg] 55 REAL(wp), PUBLIC :: tmut = 0.054_wp !: decrease of seawater meltpoint with salinity 57 56 REAL(wp), PUBLIC :: emic = 0.97_wp !: emissivity of snow or ice 58 REAL(wp), PUBLIC :: sice = 6.0_wp !: salinity of ice [psu] 59 REAL(wp), PUBLIC :: soce = 34.7_wp !: salinity of sea [psu] 60 REAL(wp), PUBLIC :: cevap = 2.5e+6_wp !: latent heat of evaporation (water) 61 REAL(wp), PUBLIC :: srgamma = 0.9_wp !: correction factor for solar radiation (Oberhuber, 1974) 57 REAL(wp), PUBLIC :: sice = 6.0_wp !: salinity of ice [psu] 58 REAL(wp), PUBLIC :: soce = 34.7_wp !: salinity of sea [psu] 62 59 REAL(wp), PUBLIC :: vkarmn = 0.4_wp !: von Karman constant 63 60 REAL(wp), PUBLIC :: stefan = 5.67e-8_wp !: Stefan-Boltzmann constant 64 61 65 #if defined key_si3 || defined key_cice 66 REAL(wp), PUBLIC :: rhoic = 917._wp !: volumic mass of sea ice [kg/m3] 67 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: thermal conductivity of fresh ice [W/m/K] 68 REAL(wp), PUBLIC :: cpic = 2067.0_wp !: specific heat of fresh ice [J/kg/K] 69 REAL(wp), PUBLIC :: lsub = 2.834e+6_wp !: pure ice latent heat of sublimation [J/kg] 70 REAL(wp), PUBLIC :: lfus = 0.334e+6_wp !: latent heat of fusion of fresh ice [J/kg] 71 REAL(wp), PUBLIC :: tmut = 0.054_wp !: decrease of seawater meltpoint with salinity 72 REAL(wp), PUBLIC :: xlsn !: = lfus*rhosn (volumetric latent heat fusion of snow) [J/m3] 73 #else 74 REAL(wp), PUBLIC :: rhoic = 900._wp !: volumic mass of sea ice [kg/m3] 75 REAL(wp), PUBLIC :: rcdic = 2.034396_wp !: conductivity of the ice [W/m/K] 76 REAL(wp), PUBLIC :: rcpic = 1.8837e+6_wp !: volumetric specific heat for ice [J/m3/K] 77 REAL(wp), PUBLIC :: cpic !: = rcpic / rhoic (specific heat for ice) [J/Kg/K] 78 REAL(wp), PUBLIC :: rcdsn = 0.22_wp !: conductivity of the snow [W/m/K] 79 REAL(wp), PUBLIC :: rcpsn = 6.9069e+5_wp !: volumetric specific heat for snow [J/m3/K] 80 REAL(wp), PUBLIC :: xlsn = 110.121e+6_wp !: volumetric latent heat fusion of snow [J/m3] 81 REAL(wp), PUBLIC :: lfus !: = xlsn / rhosn (latent heat of fusion of fresh ice) [J/Kg] 82 REAL(wp), PUBLIC :: xlic = 300.33e+6_wp !: volumetric latent heat fusion of ice [J/m3] 83 REAL(wp), PUBLIC :: xsn = 2.8e+6_wp !: volumetric latent heat of sublimation of snow [J/m3] 84 #endif 85 #if defined key_cice 86 REAL(wp), PUBLIC :: rcdsn = 0.31_wp !: thermal conductivity of snow [W/m/K] 87 #endif 88 #if defined key_si3 89 REAL(wp), PUBLIC :: r1_rhoic !: 1 / rhoic 90 REAL(wp), PUBLIC :: r1_rhosn !: 1 / rhosn 91 REAL(wp), PUBLIC :: r1_cpic !: 1 / cpic 92 #endif 62 REAL(wp), PUBLIC :: r1_rhoi !: 1 / rhoi 63 REAL(wp), PUBLIC :: r1_rhos !: 1 / rhos 64 REAL(wp), PUBLIC :: r1_rhow !: 1 / rhow 65 REAL(wp), PUBLIC :: r1_cpi !: 1 / rcpi 66 REAL(wp), PUBLIC :: r1_Lsub !: 1 / rLsub 67 REAL(wp), PUBLIC :: r1_Lfus !: 1 / rLfus 68 93 69 !!---------------------------------------------------------------------- 94 70 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 105 81 !! ** Purpose : set and print the constants 106 82 !!---------------------------------------------------------------------- 107 83 ! 108 84 IF(lwp) WRITE(numout,*) 109 85 IF(lwp) WRITE(numout,*) 'phy_cst : initialization of ocean parameters and constants' 110 86 IF(lwp) WRITE(numout,*) '~~~~~~~' 111 87 112 ! Define & print constants 113 ! ------------------------ 114 IF(lwp) WRITE(numout,*) 115 IF(lwp) WRITE(numout,*) ' Constants' 116 117 IF(lwp) WRITE(numout,*) 118 IF(lwp) WRITE(numout,*) ' mathematical constant rpi = ', rpi 88 ! !== Define derived constant ==! 119 89 120 90 rsiyea = 365.25_wp * rday * 2._wp * rpi / 6.283076_wp … … 125 95 omega = 2._wp * rpi / rsiday 126 96 #endif 127 IF(lwp) WRITE(numout,*) 128 IF(lwp) WRITE(numout,*) ' day rday = ', rday, ' s' 129 IF(lwp) WRITE(numout,*) ' sideral year rsiyea = ', rsiyea, ' s' 130 IF(lwp) WRITE(numout,*) ' sideral day rsiday = ', rsiday, ' s' 131 IF(lwp) WRITE(numout,*) ' omega omega = ', omega, ' s^-1' 132 IF(lwp) WRITE(numout,*) 133 IF(lwp) WRITE(numout,*) ' nb of months per year raamo = ', raamo, ' months' 134 IF(lwp) WRITE(numout,*) ' nb of hours per day rjjhh = ', rjjhh, ' hours' 135 IF(lwp) WRITE(numout,*) ' nb of minutes per hour rhhmm = ', rhhmm, ' mn' 136 IF(lwp) WRITE(numout,*) ' nb of seconds per minute rmmss = ', rmmss, ' s' 137 IF(lwp) WRITE(numout,*) 138 IF(lwp) WRITE(numout,*) ' earth radius ra = ', ra, ' m' 139 IF(lwp) WRITE(numout,*) ' gravity grav = ', grav , ' m/s^2' 140 IF(lwp) WRITE(numout,*) 141 IF(lwp) WRITE(numout,*) ' triple point of temperature rtt = ', rtt , ' K' 142 IF(lwp) WRITE(numout,*) ' freezing point of water rt0 = ', rt0 , ' K' 143 IF(lwp) WRITE(numout,*) ' melting point of snow rt0_snow = ', rt0_snow, ' K' 144 IF(lwp) WRITE(numout,*) ' melting point of ice rt0_ice = ', rt0_ice , ' K' 145 IF(lwp) WRITE(numout,*) 146 IF(lwp) WRITE(numout,*) ' reference density and heat capacity now defined in eosbn2.f90' 147 148 #if defined key_si3 || defined key_cice 149 xlsn = lfus * rhosn ! volumetric latent heat fusion of snow [J/m3] 150 #else 151 cpic = rcpic / rhoic ! specific heat for ice [J/Kg/K] 152 lfus = xlsn / rhosn ! latent heat of fusion of fresh ice 153 #endif 154 #if defined key_si3 155 r1_rhoic = 1._wp / rhoic 156 r1_rhosn = 1._wp / rhosn 157 r1_cpic = 1._wp / cpic 158 #endif 159 IF(lwp) THEN 97 98 r1_rhoi = 1._wp / rhoi 99 r1_rhos = 1._wp / rhos 100 r1_cpi = 1._wp / rcpi 101 r1_Lsub = 1._wp / rLsub 102 r1_Lfus = 1._wp / rLfus 103 104 IF(lwp) THEN !== print constants ==! 160 105 WRITE(numout,*) 161 #if defined key_cice 162 WRITE(numout,*) ' thermal conductivity of the snow = ', rcdsn , ' J/s/m/K' 163 #endif 164 WRITE(numout,*) ' thermal conductivity of pure ice = ', rcdic , ' J/s/m/K' 165 WRITE(numout,*) ' fresh ice specific heat = ', cpic , ' J/kg/K' 166 WRITE(numout,*) ' latent heat of fusion of fresh ice / snow = ', lfus , ' J/kg' 167 #if defined key_si3 || defined key_cice 168 WRITE(numout,*) ' latent heat of subl. of fresh ice / snow = ', lsub , ' J/kg' 169 #else 170 WRITE(numout,*) ' density times specific heat for snow = ', rcpsn , ' J/m^3/K' 171 WRITE(numout,*) ' density times specific heat for ice = ', rcpic , ' J/m^3/K' 172 WRITE(numout,*) ' volumetric latent heat fusion of sea ice = ', xlic , ' J/m' 173 WRITE(numout,*) ' latent heat of sublimation of snow = ', xsn , ' J/kg' 174 #endif 175 WRITE(numout,*) ' volumetric latent heat fusion of snow = ', xlsn , ' J/m^3' 176 WRITE(numout,*) ' density of sea ice = ', rhoic , ' kg/m^3' 177 WRITE(numout,*) ' density of snow = ', rhosn , ' kg/m^3' 178 WRITE(numout,*) ' density of freshwater (in melt ponds) = ', rhofw , ' kg/m^3' 179 WRITE(numout,*) ' emissivity of snow or ice = ', emic 180 WRITE(numout,*) ' salinity of ice = ', sice , ' psu' 181 WRITE(numout,*) ' salinity of sea = ', soce , ' psu' 182 WRITE(numout,*) ' latent heat of evaporation (water) = ', cevap , ' J/m^3' 183 WRITE(numout,*) ' correction factor for solar radiation = ', srgamma 184 WRITE(numout,*) ' von Karman constant = ', vkarmn 185 WRITE(numout,*) ' Stefan-Boltzmann constant = ', stefan , ' J/s/m^2/K^4' 106 WRITE(numout,*) ' Constants' 186 107 WRITE(numout,*) 187 WRITE(numout,*) ' conversion: degre ==> radian rad = ', rad 108 WRITE(numout,*) ' mathematical constant rpi = ', rpi 109 WRITE(numout,*) ' conversion: degre ==> radian rad = ', rad 188 110 WRITE(numout,*) 189 WRITE(numout,*) ' smallest real computer value rsmall = ', rsmall 111 WRITE(numout,*) ' day in seconds rday = ', rday , ' s' 112 WRITE(numout,*) ' sideral year rsiyea = ', rsiyea, ' s' 113 WRITE(numout,*) ' sideral day rsiday = ', rsiday, ' s' 114 WRITE(numout,*) ' omega = 2 pi / rsiday omega = ', omega , ' s^-1' 115 WRITE(numout,*) ' earth radius ra = ', ra , ' m' 116 WRITE(numout,*) ' gravity grav = ', grav , ' m/s^2' 117 WRITE(numout,*) 118 WRITE(numout,*) ' nb of months per year raamo = ', raamo, ' months' 119 WRITE(numout,*) ' nb of hours per day rjjhh = ', rjjhh, ' hours' 120 WRITE(numout,*) ' nb of minutes per hour rhhmm = ', rhhmm, ' mn' 121 WRITE(numout,*) ' nb of seconds per minute rmmss = ', rmmss, ' s' 122 WRITE(numout,*) 123 WRITE(numout,*) ' reference ocean density and heat capacity now defined in eosbn2.f90' 124 WRITE(numout,*) 125 WRITE(numout,*) ' freezing point of freshwater rt0 = ', rt0 , ' K' 126 WRITE(numout,*) ' sea ice density rhoi = ', rhoi , ' kg/m^3' 127 WRITE(numout,*) ' snow density rhos = ', rhos , ' kg/m^3' 128 WRITE(numout,*) ' freshwater density (in melt ponds) rhow = ', rhow , ' kg/m^3' 129 WRITE(numout,*) ' thermal conductivity of pure ice rcnd_i = ', rcnd_i, ' J/s/m/K' 130 WRITE(numout,*) ' fresh ice specific heat rcpi = ', rcpi , ' J/kg/K' 131 WRITE(numout,*) ' latent heat of fusion of fresh ice / snow rLfus = ', rLfus , ' J/kg' 132 WRITE(numout,*) ' latent heat of subl. of fresh ice / snow rLsub = ', rLsub , ' J/kg' 133 WRITE(numout,*) ' emissivity of snow or ice emic = ', emic 134 WRITE(numout,*) ' salinity of ice sice = ', sice , ' psu' 135 WRITE(numout,*) ' salinity of sea soce = ', soce , ' psu' 136 WRITE(numout,*) ' von Karman constant vkarmn = ', vkarmn 137 WRITE(numout,*) ' Stefan-Boltzmann constant stefan = ', stefan, ' J/s/m^2/K^4' 138 WRITE(numout,*) 139 WRITE(numout,*) 140 WRITE(numout,*) ' smallest real computer value rsmall = ', rsmall 190 141 ENDIF 191 142 ! 192 143 END SUBROUTINE phy_cst 193 144 -
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/restart.F90
r9838 r9939 8 8 !! 2.0 ! 2006-07 (S. Masson) use IOM for restart 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 10 !! - - ! 2010-10 (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D) 11 !! 3.7 ! 2014-01 (G. Madec) suppression of curl and hdiv from the restart 12 !! - ! 2014-12 (G. Madec) remove KPP scheme 13 !!---------------------------------------------------------------------- 14 15 !!---------------------------------------------------------------------- 16 !! rst_opn : open the ocean restart file 17 !! rst_write : write the ocean restart file 18 !! rst_read : read the ocean restart file 19 !!---------------------------------------------------------------------- 20 USE oce ! ocean dynamics and tracers 21 USE dom_oce ! ocean space and time domain 22 USE sbc_ice ! only lk_si3 23 USE phycst ! physical constants 24 USE eosbn2 ! equation of state (eos bn2 routine) 25 USE trdmxl_oce ! ocean active mixed layer tracers trends variables 10 !! - - ! 2010-10 (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D) 11 !! 3.7 ! 2014-01 (G. Madec) suppression of curl and hdiv from the restart 12 !! - ! 2014-12 (G. Madec) remove KPP scheme 13 !! 4.0 ! 2018-06 (G. Madec) introduce l_1st_euler 14 !!---------------------------------------------------------------------- 15 16 !!---------------------------------------------------------------------- 17 !! rst_opn : open the ocean restart file in write mode 18 !! rst_write : write the ocean restart file 19 !! rst_read_open : open the ocean restart file in read mode 20 !! rst_read : read the ocean restart file 21 !!---------------------------------------------------------------------- 22 USE oce ! ocean dynamics and tracers 23 USE dom_oce ! ocean space and time domain 24 USE sbc_ice ! only lk_si3 25 USE phycst ! physical constants 26 USE eosbn2 ! equation of state (eos bn2 routine) 27 USE trdmxl_oce ! ocean active mixed layer tracers trends variables 28 USE diurnal_bulk ! 26 29 ! 27 USE in_out_manager ! I/O manager 28 USE iom ! I/O module 29 USE diurnal_bulk 30 USE in_out_manager ! I/O manager 31 USE iom ! I/O module 30 32 31 33 IMPLICIT NONE … … 34 36 PUBLIC rst_opn ! routine called by step module 35 37 PUBLIC rst_write ! routine called by step module 38 PUBLIC rst_read_open ! routine called in rst_read and (possibly) in dom_vvl_init 36 39 PUBLIC rst_read ! routine called by istate module 37 PUBLIC rst_read_open ! routine called in rst_read and (possibly) in dom_vvl_init38 40 39 41 !! * Substitutions … … 144 146 INTEGER, INTENT(in) :: kt ! ocean time-step 145 147 !!---------------------------------------------------------------------- 146 IF(lwxios) CALL iom_swap( cwxios_context ) 147 CALL iom_rstput( kt, nitrst, numrow, 'rdt' , rdt , ldxios = lwxios) ! dynamics time step 148 149 IF ( .NOT. ln_diurnal_only ) THEN 150 CALL iom_rstput( kt, nitrst, numrow, 'ub' , ub, ldxios = lwxios ) ! before fields 151 CALL iom_rstput( kt, nitrst, numrow, 'vb' , vb, ldxios = lwxios ) 152 CALL iom_rstput( kt, nitrst, numrow, 'tb' , tsb(:,:,:,jp_tem), ldxios = lwxios ) 153 CALL iom_rstput( kt, nitrst, numrow, 'sb' , tsb(:,:,:,jp_sal), ldxios = lwxios ) 154 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb, ldxios = lwxios ) 155 ! 156 CALL iom_rstput( kt, nitrst, numrow, 'un' , un, ldxios = lwxios ) ! now fields 157 CALL iom_rstput( kt, nitrst, numrow, 'vn' , vn, ldxios = lwxios ) 158 CALL iom_rstput( kt, nitrst, numrow, 'tn' , tsn(:,:,:,jp_tem), ldxios = lwxios ) 159 CALL iom_rstput( kt, nitrst, numrow, 'sn' , tsn(:,:,:,jp_sal), ldxios = lwxios ) 160 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn, ldxios = lwxios ) 161 CALL iom_rstput( kt, nitrst, numrow, 'rhop' , rhop, ldxios = lwxios ) 162 ! extra variable needed for the ice sheet coupling 163 IF ( ln_iscpl ) THEN 164 CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask, ldxios = lwxios ) ! need to extrapolate T/S 165 CALL iom_rstput( kt, nitrst, numrow, 'umask' , umask, ldxios = lwxios ) ! need to correct barotropic velocity 166 CALL iom_rstput( kt, nitrst, numrow, 'vmask' , vmask, ldxios = lwxios ) ! need to correct barotropic velocity 167 CALL iom_rstput( kt, nitrst, numrow, 'smask' , ssmask, ldxios = lwxios) ! need to correct barotropic velocity 168 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios ) ! need to compute temperature correction 169 CALL iom_rstput( kt, nitrst, numrow, 'e3u_n', e3u_n(:,:,:), ldxios = lwxios ) ! need to compute bt conservation 170 CALL iom_rstput( kt, nitrst, numrow, 'e3v_n', e3v_n(:,:,:), ldxios = lwxios ) ! need to compute bt conservation 171 CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n(:,:,:), ldxios = lwxios ) ! need to compute extrapolation if vvl 172 END IF 173 ENDIF 174 175 IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst, ldxios = lwxios ) 176 IF(lwxios) CALL iom_swap( cxios_context ) 148 IF( lwxios ) CALL iom_swap( cwxios_context ) 149 150 CALL iom_rstput( kt, nitrst, numrow, 'rdt' , rn_Dt , ldxios = lwxios ) ! dynamics time step 151 ! 152 IF( .NOT. ln_diurnal_only ) THEN 153 CALL iom_rstput( kt, nitrst, numrow, 'ub' , ub , ldxios = lwxios ) ! before fields 154 CALL iom_rstput( kt, nitrst, numrow, 'vb' , vb , ldxios = lwxios ) 155 CALL iom_rstput( kt, nitrst, numrow, 'tb' , tsb(:,:,:,jp_tem), ldxios = lwxios ) 156 CALL iom_rstput( kt, nitrst, numrow, 'sb' , tsb(:,:,:,jp_sal), ldxios = lwxios ) 157 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb , ldxios = lwxios ) 158 ! 159 CALL iom_rstput( kt, nitrst, numrow, 'un' , un , ldxios = lwxios ) ! now fields 160 CALL iom_rstput( kt, nitrst, numrow, 'vn' , vn , ldxios = lwxios ) 161 CALL iom_rstput( kt, nitrst, numrow, 'tn' , tsn(:,:,:,jp_tem), ldxios = lwxios ) 162 CALL iom_rstput( kt, nitrst, numrow, 'sn' , tsn(:,:,:,jp_sal), ldxios = lwxios ) 163 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn , ldxios = lwxios ) 164 CALL iom_rstput( kt, nitrst, numrow, 'rhop' , rhop , ldxios = lwxios ) 165 ! 166 IF( ln_iscpl ) THEN ! extra variable needed for the ice sheet coupling 167 CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask , ldxios = lwxios ) ! need to extrapolate T/S 168 CALL iom_rstput( kt, nitrst, numrow, 'umask' , umask , ldxios = lwxios ) ! need to correct barotropic velocity 169 CALL iom_rstput( kt, nitrst, numrow, 'vmask' , vmask , ldxios = lwxios ) ! need to correct barotropic velocity 170 CALL iom_rstput( kt, nitrst, numrow, 'smask' , ssmask , ldxios = lwxios ) ! need to correct barotropic velocity 171 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n' , e3t_n , ldxios = lwxios ) ! need to compute temperature correction 172 CALL iom_rstput( kt, nitrst, numrow, 'e3u_n' , e3u_n , ldxios = lwxios ) ! need to compute bt conservation 173 CALL iom_rstput( kt, nitrst, numrow, 'e3v_n' , e3v_n , ldxios = lwxios ) ! need to compute bt conservation 174 CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n, ldxios = lwxios ) ! need to compute extrapolation if vvl 175 ENDIF 176 ENDIF 177 ! 178 IF( ln_diurnal ) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst, ldxios = lwxios ) 179 IF( lwxios ) CALL iom_swap( cxios_context ) 177 180 IF( kt == nitrst ) THEN 178 IF(.NOT.lwxios) THEN 179 CALL iom_close( numrow ) ! close the restart file (only at last time step) 180 ELSE 181 CALL iom_context_finalize( cwxios_context ) 181 IF( lwxios ) THEN ; CALL iom_context_finalize( cwxios_context ) 182 ELSE ; CALL iom_close( numrow ) ! close the restart file (only at last time step) 182 183 ENDIF 183 184 !!gm IF( .NOT. lk_trdmld ) lrst_oce = .FALSE. 184 185 !!gm not sure what to do here ===>>> ask to Sebastian 185 186 lrst_oce = .FALSE. 186 187 188 189 187 IF( ln_rst_list ) THEN 188 nrst_lst = MIN(nrst_lst + 1, SIZE(nstocklist,1)) 189 nitrst = nstocklist( nrst_lst ) 190 ENDIF 190 191 ENDIF 191 192 ! … … 202 203 !! the file has already been opened 203 204 !!---------------------------------------------------------------------- 204 INTEGER 205 LOGICAL 206 CHARACTER(lc) 205 INTEGER :: jlibalt = jprstlib 206 LOGICAL :: llok 207 CHARACTER(lc) :: clpath ! full path to ocean output restart file 207 208 !!---------------------------------------------------------------------- 208 209 ! … … 238 239 ENDIF 239 240 ENDIF 240 241 ! 241 242 END SUBROUTINE rst_read_open 242 243 … … 254 255 REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d 255 256 !!---------------------------------------------------------------------- 256 257 ! 257 258 CALL rst_read_open ! open restart for reading (if not already opened) 258 259 … … 260 261 IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 ) THEN 261 262 CALL iom_get( numror, 'rdt', zrdt, ldxios = lrxios ) 262 IF( zrdt /= rdt ) neuler = 0 263 IF( zrdt /= rn_Dt ) THEN 264 IF(lwp) WRITE( numout,*) 265 IF(lwp) WRITE( numout,*) 'rst_read: rdt not equal to the read one' 266 IF(lwp) WRITE( numout,*) 267 IF(lwp) WRITE( numout,*) ' ==>>> forced euler first time-step' 268 l_1st_euler = .TRUE. 269 ENDIF 263 270 ENDIF 264 271 … … 266 273 IF( ln_diurnal ) CALL iom_get( numror, jpdom_autoglo, 'Dsst' , x_dsst, ldxios = lrxios ) 267 274 IF ( ln_diurnal_only ) THEN 268 IF(lwp) WRITE( numout, * ) & 269 & "rst_read:- ln_diurnal_only set, setting rhop=rau0" 270 rhop = rau0 275 IF(lwp) WRITE( numout,*) 'rst_read: ln_diurnal_only set, setting rhop=rho0' 276 rhop = rho0 271 277 CALL iom_get( numror, jpdom_autoglo, 'tn' , w3d, ldxios = lrxios ) 272 278 tsn(:,:,1,jp_tem) = w3d(:,:,1) … … 274 280 ENDIF 275 281 276 IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 ) THEN282 IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 .AND. .NOT.l_1st_euler ) THEN 277 283 CALL iom_get( numror, jpdom_autoglo, 'ub' , ub, ldxios = lrxios ) ! before fields 278 284 CALL iom_get( numror, jpdom_autoglo, 'vb' , vb, ldxios = lrxios ) … … 281 287 CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb, ldxios = lrxios ) 282 288 ELSE 283 neuler = 0289 l_1st_euler = .TRUE. ! before field not found, forced euler 1st time-step 284 290 ENDIF 285 291 ! … … 295 301 ENDIF 296 302 ! 297 IF( neuler == 0 ) THEN ! Euler restart (neuler=0) 298 tsb (:,:,:,:) = tsn (:,:,:,:) ! all before fields set to now values 299 ub (:,:,:) = un (:,:,:) 300 vb (:,:,:) = vn (:,:,:) 301 sshb (:,:) = sshn (:,:) 302 ! 303 IF( .NOT.ln_linssh ) THEN 304 DO jk = 1, jpk 305 e3t_b(:,:,jk) = e3t_n(:,:,jk) 306 END DO 307 ENDIF 308 ! 303 IF( l_1st_euler ) THEN ! Euler restart 304 tsb (:,:,:,:) = tsn (:,:,:,:) ! all before fields set to now values 305 ub (:,:,:) = un (:,:,:) 306 vb (:,:,:) = vn (:,:,:) 307 sshb(:,:) = sshn(:,:) 308 IF( .NOT.ln_linssh ) e3t_b(:,:,:) = e3t_n(:,:,:) 309 309 ENDIF 310 310 !
Note: See TracChangeset
for help on using the changeset viewer.
