Changeset 9923
- Timestamp:
- 2018-07-11T10:24:17+02:00 (3 years ago)
- Location:
- NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src
- Files:
-
- 135 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/ice.F90
r9604 r9923 188 188 ! !!** some other parameters 189 189 INTEGER , PUBLIC :: kt_ice !: iteration number 190 REAL(wp), PUBLIC :: r dt_ice !: ice time step191 REAL(wp), PUBLIC :: r1_ rdtice !: = 1. / rdt_ice190 REAL(wp), PUBLIC :: rDt_ice !: ice time step 191 REAL(wp), PUBLIC :: r1_Dt_ice !: = 1. / rdt_ice 192 192 REAL(wp), PUBLIC :: r1_nlay_i !: 1 / nlay_i 193 193 REAL(wp), PUBLIC :: r1_nlay_s !: 1 / nlay_s -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icecor.F90
r9604 r9923 86 86 IF ( nn_icesal == 2 ) THEN ! salinity must stay in bounds [Simin,Simax] ! 87 87 ! !----------------------------------------------------- 88 zzc = rhoic * r1_ rdtice88 zzc = rhoic * r1_Dt_ice 89 89 DO jl = 1, jpl 90 90 DO jj = 1, jpj … … 137 137 ! ! heat content variation (W.m-2) 138 138 diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) & 139 & + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) ) * r1_ rdtice139 & + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) ) * r1_Dt_ice 140 140 ! ! salt, volume 141 diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_ rdtice142 diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_ rdtice143 diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_ rdtice141 diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_Dt_ice 142 diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_Dt_ice 143 diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_Dt_ice 144 144 END DO 145 145 END DO 146 146 ! ! concentration tendency (dynamics) 147 zafx (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_ rdtice147 zafx (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 148 148 afx_tot(:,:) = zafx(:,:) 149 149 IF( iom_use('afxdyn') ) CALL iom_put( 'afxdyn' , zafx(:,:) ) … … 158 158 ! ! heat content variation (W.m-2) 159 159 diag_heat(ji,jj) = diag_heat(ji,jj) - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) & 160 & + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) ) * r1_ rdtice160 & + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) ) * r1_Dt_ice 161 161 ! ! salt, volume 162 diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_ rdtice163 diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_ rdtice164 diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_ rdtice162 diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_Dt_ice 163 diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_Dt_ice 164 diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_Dt_ice 165 165 END DO 166 166 END DO 167 167 ! ! concentration tendency (total + thermo) 168 zafx (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_ rdtice168 zafx (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 169 169 afx_tot(:,:) = afx_tot(:,:) + zafx(:,:) 170 170 IF( iom_use('afxthd') ) CALL iom_put( 'afxthd' , zafx(:,:) ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icectl.F90
r9604 r9923 121 121 ! outputs 122 122 zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv & 123 & - pdiag_v ) * r1_ rdtice - zfv ) * rday123 & - pdiag_v ) * r1_Dt_ice - zfv ) * rday 124 124 125 125 zs = ( ( glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t ) * zconv & 126 & - pdiag_s ) * r1_ rdtice + zfs ) * rday126 & - pdiag_s ) * r1_Dt_ice + zfs ) * rday 127 127 128 128 zt = ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) & 129 129 & + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv & 130 & - pdiag_t ) * r1_ rdtice + zft130 & - pdiag_t ) * r1_Dt_ice + zft 131 131 132 132 ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative … … 580 580 WRITE(numout,*) ' hfx_res : ', hfx_res(ji,jj) 581 581 WRITE(numout,*) ' fhtur : ', fhtur(ji,jj) 582 WRITE(numout,*) ' qlead : ', qlead(ji,jj) * r1_ rdtice582 WRITE(numout,*) ' qlead : ', qlead(ji,jj) * r1_Dt_ice 583 583 WRITE(numout,*) 584 584 WRITE(numout,*) ' - Salt fluxes at bottom interface ***' -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedia.F90
r9604 r9923 95 95 ! 2 - Trends due to forcing ! 96 96 ! ---------------------------! 97 z_frc_volbot = r1_r au0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater flux ice/snow-ocean98 z_frc_voltop = r1_r au0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater flux ice/snow-atm99 z_frc_sal = r1_r au0 * glob_sum( - sfx(:,:) * e1e2t(:,:) ) * 1.e-9 ! salt fluxes ice/snow-ocean97 z_frc_volbot = r1_rho0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater flux ice/snow-ocean 98 z_frc_voltop = r1_rho0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater flux ice/snow-atm 99 z_frc_sal = r1_rho0 * glob_sum( - sfx(:,:) * e1e2t(:,:) ) * 1.e-9 ! salt fluxes ice/snow-ocean 100 100 z_frc_tembot = glob_sum( hfx_out(:,:) * e1e2t(:,:) ) * 1.e-20 ! heat on top of ocean (and below ice) 101 101 z_frc_temtop = glob_sum( hfx_in (:,:) * e1e2t(:,:) ) * 1.e-20 ! heat on top of ice-coean … … 110 110 ! 3 - Content variations ! 111 111 ! ----------------------- ! 112 zdiff_vol = r1_r au0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3)113 zdiff_sal = r1_r au0 * glob_sum( ( rhoic* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss)112 zdiff_vol = r1_rho0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3) 113 zdiff_sal = r1_rho0 * glob_sum( ( rhoic* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss) 114 114 zdiff_tem = glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J) 115 115 ! + SUM( qevap_ice * a_i_b, dim=3 ) !! clem: I think this term should not be there (but needs a check) … … 125 125 ! 5 - Diagnostics writing ! 126 126 ! ----------------------- ! 127 !!gm I don't understand the division by the ocean surface (i.e. glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt )128 !! and its multiplication bu kt ! is it really what we want ? what is this quantity ?129 !! IF it is really what we want, compute it at kt=nit000, not 3 time by time-step !130 !! kt*rdt : you mean rdtice ?131 !!gm132 127 ! 133 128 IF( iom_use('ibgvolume') ) CALL iom_put( 'ibgvolume' , zdiff_vol ) ! ice/snow volume drift (km3 equivalent ocean water) … … 135 130 IF( iom_use('ibgheatco') ) CALL iom_put( 'ibgheatco' , zdiff_tem ) ! ice/snow heat content drift (1.e20 J) 136 131 IF( iom_use('ibgheatfx') ) CALL iom_put( 'ibgheatfx' , & ! ice/snow heat flux drift (W/m2) 137 & zdiff_tem /glob_sum( e1e2t(:,:) * 1.e-20 * kt*r dt ) )132 & zdiff_tem /glob_sum( e1e2t(:,:) * 1.e-20 * kt*rn_Dt ) ) 138 133 139 134 IF( iom_use('ibgfrcvoltop') ) CALL iom_put( 'ibgfrcvoltop' , frc_voltop ) ! vol forcing ice/snw-atm (km3 equivalent ocean water) … … 143 138 IF( iom_use('ibgfrctembot') ) CALL iom_put( 'ibgfrctembot' , frc_tembot ) ! heat on top of ocean(below ice) (1.e20 J) 144 139 IF( iom_use('ibgfrchfxtop') ) CALL iom_put( 'ibgfrchfxtop' , & ! heat on top of ice/snw/ocean (W/m2) 145 & frc_temtop / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*r dt )140 & frc_temtop / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rn_Dt ) 146 141 IF( iom_use('ibgfrchfxbot') ) CALL iom_put( 'ibgfrchfxbot' , & ! heat on top of ocean(below ice) (W/m2) 147 & frc_tembot / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*r dt )142 & frc_tembot / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rn_Dt ) 148 143 149 144 IF( iom_use('ibgvol_tot' ) ) CALL iom_put( 'ibgvol_tot' , zbg_ivol ) ! ice volume (km3) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_adv.F90
r9604 r9923 98 98 ! diagnostics 99 99 !------------ 100 diag_trp_ei(:,:) = SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_ rdtice101 diag_trp_es(:,:) = SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_ rdtice102 diag_trp_sv(:,:) = SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_ rdtice103 diag_trp_vi(:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_ rdtice104 diag_trp_vs(:,:) = SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_ rdtice100 diag_trp_ei(:,:) = SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice 101 diag_trp_es(:,:) = SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 102 diag_trp_sv(:,:) = SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice 103 diag_trp_vi(:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice 104 diag_trp_vs(:,:) = SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice 105 105 IF( iom_use('icemtrp') ) CALL iom_put( "icemtrp" , diag_trp_vi * rhoic ) ! ice mass transport 106 106 IF( iom_use('snwmtrp') ) CALL iom_put( "snwmtrp" , diag_trp_vs * rhosn ) ! snw mass transport -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rdgrft.F90
r9604 r9923 189 189 ! divergence given by the advection scheme 190 190 ! (which may not be equal to divu as computed from the velocity field) 191 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_ rdtice191 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_Dt_ice 192 192 ! 193 193 IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) ! make sure the closing rate is large enough … … 255 255 ELSE 256 256 iterate_ridging = 1 257 zdivu_adv (ji) = zfac * r1_ rdtice257 zdivu_adv (ji) = zfac * r1_Dt_ice 258 258 closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) ) 259 259 opning (ji) = MAX( 0._wp, zdivu_adv(ji) ) … … 460 460 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 461 461 IF( zfac > pa_i(ji,jl) ) THEN 462 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_ rdtice462 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice 463 463 ENDIF 464 464 END DO … … 472 472 zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice 473 473 IF( zfac < 0._wp ) THEN ! would lead to negative ato_i 474 opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_ rdtice474 opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice 475 475 ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum 476 opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_ rdtice476 opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice 477 477 ENDIF 478 478 END DO … … 570 570 571 571 ! Ice-ocean exchanges associated with ice porosity 572 wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids573 sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_ rdtice574 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice! > 0 [W.m-2]572 wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_Dt_ice ! increase in ice volume due to seawater frozen in voids 573 sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_Dt_ice 574 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_Dt_ice ! > 0 [W.m-2] 575 575 576 576 ! Put the snow lost by ridging into the ocean 577 577 ! Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 578 578 wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! fresh water source for ocean 579 & + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_ rdtice579 & + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice 580 580 581 581 ! Put the melt pond water into the ocean … … 584 584 !IF ( ln_pnd_fwb ) THEN 585 585 ! wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg ) & ! fresh water source for ocean 586 ! & + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_ rdtice586 ! & + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_Dt_ice 587 587 !ENDIF 588 588 … … 590 590 IF( nn_icesal /= 2 ) THEN 591 591 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i 592 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_ rdtice & ! put back sss_m into the ocean593 & - s_i_1d(ji) * vsw * rhoic * r1_ rdtice ! and get s_i from the ocean592 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_Dt_ice & ! put back sss_m into the ocean 593 & - s_i_1d(ji) * vsw * rhoic * r1_Dt_ice ! and get s_i from the ocean 594 594 ENDIF 595 595 … … 621 621 ! Put the snow lost by ridging into the ocean 622 622 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2) 623 & - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_ rdtice623 & - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice 624 624 ! 625 625 ! Remove energy of new ridge to each category jl1 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rhg_evp.F90
r9660 r9923 114 114 INTEGER :: jter ! local integers 115 115 ! 116 REAL(wp) :: zrhoco ! r au0 * rn_cio116 REAL(wp) :: zrhoco ! rho0 * rn_cio 117 117 REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling 118 118 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity … … 221 221 ! 1) define some variables and initialize arrays 222 222 !------------------------------------------------------------------------------! 223 zrhoco = r au0 * rn_cio223 zrhoco = rho0 * rn_cio 224 224 225 225 ! ecc2: square of yield ellipse eccenticrity … … 271 271 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 272 272 ! 273 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_r au0273 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rho0 274 274 ! 275 275 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceistate.F90
r9656 r9923 295 295 ! In case snow load is in excess that would lead to transformation from snow to ice 296 296 ! Then, transfer the snow excess into the ice (different from icethd_dh) 297 zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - r au0 ) * h_i(ji,jj,jl) ) * r1_rau0 )297 zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 ) 298 298 ! recompute h_i, h_s avoiding out of bounds values 299 299 h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) … … 415 415 IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area 416 416 ! 417 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_r au0418 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_r au0417 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rho0 418 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rho0 419 419 ! 420 420 IF( .NOT.ln_linssh ) THEN -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icestp.F90
r9725 r9923 341 341 IF( ln_bdy .AND. ln_icediachk ) CALL ctl_warn('par_init: online conservation check does not work with BDY') 342 342 ! 343 rdt_ice = REAL(nn_fsbc) * r dt!--- sea-ice timestep and its inverse344 r1_ rdtice = 1._wp / rdt_ice343 rdt_ice = REAL(nn_fsbc) * rn_Dt !--- sea-ice timestep and its inverse 344 r1_Dt_ice = 1._wp / rdt_ice 345 345 IF(lwp) WRITE(numout,*) 346 IF(lwp) WRITE(numout,*) ' ice timestep rdt_ice = nn_fsbc*r dt = ', rdt_ice346 IF(lwp) WRITE(numout,*) ' ice timestep rdt_ice = nn_fsbc*rn_Dt = ', rdt_ice, ' [s]' 347 347 ! 348 348 r1_nlay_i = 1._wp / REAL( nlay_i, wp ) !--- inverse of nlay_i and nlay_s -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd.F90
r9750 r9923 120 120 DO jj = 2, jpjm1 121 121 DO ji = fs_2, fs_jpim1 122 zfric(ji,jj) = r1_r au0 * SQRT( 0.5_wp * &122 zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp * & 123 123 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 124 124 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) … … 150 150 ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 151 151 ! includes supercooling potential energy (>0) or "above-freezing" energy (<0) 152 zqfr = tmask(ji,jj,1) * r au0 *rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )152 zqfr = tmask(ji,jj,1) * rho0_rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 153 153 154 154 ! --- Above-freezing sensible heat content (J/m2 grid) 155 zqfr_neg = tmask(ji,jj,1) * r au0 *rcp * e3t_m(ji,jj) * MIN( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0._wp )155 zqfr_neg = tmask(ji,jj,1) * rho0_rcp * e3t_m(ji,jj) * MIN( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0._wp ) 156 156 157 157 ! --- Sensible ocean-to-ice heat flux (W/m2) 158 158 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 159 fhtur(ji,jj) = rswitch * r au0 *rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2160 161 fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_neg * r1_ rdtice / MAX( at_i(ji,jj), epsi10 ) )159 fhtur(ji,jj) = rswitch * rho0_rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 160 161 fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 162 162 ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach 163 163 ! the freezing point, so that we do not have SST < T_freeze … … 169 169 ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 170 170 IF( zqld > 0._wp ) THEN 171 fhld (ji,jj) = rswitch * zqld * r1_ rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90171 fhld (ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 172 172 qlead(ji,jj) = 0._wp 173 173 ELSE … … 197 197 ! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar 198 198 hfx_out(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean 199 & - qlead(:,:) * r1_ rdtice & ! heat flux taken from the ocean where there is open water ice formation199 & - qlead(:,:) * r1_Dt_ice & ! heat flux taken from the ocean where there is open water ice formation 200 200 & - at_i (:,:) * fhtur(:,:) & ! heat flux taken by turbulence 201 201 & - at_i (:,:) * fhld(:,:) ! heat flux taken during bottom growth/melt -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_da.F90
r9604 r9923 137 137 138 138 ! Contribution to salt flux 139 sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoic * h_i_1d(ji) * zda * s_i_1d(ji) * r1_ rdtice139 sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoic * h_i_1d(ji) * zda * s_i_1d(ji) * r1_Dt_ice 140 140 141 141 ! Contribution to heat flux into the ocean [W.m-2], (<0) 142 hfx_thd_1d(ji) = hfx_thd_1d(ji) - zda * r1_ rdtice * ( h_i_1d(ji) * r1_nlay_i * SUM( e_i_1d(ji,1:nlay_i) ) &142 hfx_thd_1d(ji) = hfx_thd_1d(ji) - zda * r1_Dt_ice * ( h_i_1d(ji) * r1_nlay_i * SUM( e_i_1d(ji,1:nlay_i) ) & 143 143 + h_s_1d(ji) * r1_nlay_s * SUM( e_s_1d(ji,1:nlay_s) ) ) 144 144 145 145 ! Contribution to mass flux 146 wfx_lam_1d(ji) = wfx_lam_1d(ji) + zda * r1_ rdtice * ( rhoic * h_i_1d(ji) + rhosn * h_s_1d(ji) )146 wfx_lam_1d(ji) = wfx_lam_1d(ji) + zda * r1_Dt_ice * ( rhoic * h_i_1d(ji) + rhosn * h_s_1d(ji) ) 147 147 148 148 ! new concentration -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_dh.F90
r9767 r9923 76 76 REAL(wp) :: zgrr ! bottom growth rate 77 77 REAL(wp) :: zt_i_new ! bottom formation temperature 78 REAL(wp) :: z1_rho ! 1/(rhosn+r au0-rhoic)78 REAL(wp) :: z1_rho ! 1/(rhosn+rho0-rhoic) 79 79 80 80 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean … … 181 181 DO ji = 1, npti 182 182 IF( t_s_1d(ji,jk) > rt0 ) THEN 183 hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_ rdtice ! heat flux to the ocean [W.m-2], < 0184 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * zh_s(ji,jk) * a_i_1d(ji) * r1_ rdtice ! mass flux183 hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0 184 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux 185 185 ! updates 186 186 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk) … … 205 205 zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3) 206 206 ! 207 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_ rdtice ! heat flux from snow precip (>0, W.m-2)208 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_ rdtice ! mass flux, <0207 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_Dt_ice ! heat flux from snow precip (>0, W.m-2) 208 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice ! mass flux, <0 209 209 210 210 ! --- melt of falling snow --- … … 212 212 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change 213 213 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 214 hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_ rdtice ! heat used to melt snow (W.m-2, >0)215 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_ rdtice ! snow melting only = water into the ocean (then without snow precip), >0214 hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_Dt_ice ! heat used to melt snow (W.m-2, >0) 215 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! snow melting only = water into the ocean (then without snow precip), >0 216 216 217 217 ! updates available heat + precipitations after melting … … 252 252 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 253 253 254 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_ rdtice ! heat used to melt snow(W.m-2, >0)255 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_ rdtice ! snow melting only = water into the ocean (then without snow precip)254 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0) 255 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! snow melting only = water into the ocean (then without snow precip) 256 256 257 257 ! updates available heat + thickness … … 279 279 hfx_sub_1d (ji) = hfx_sub_1d(ji) + & ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 280 280 & ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) & 281 & * a_i_1d(ji) * r1_ rdtice282 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_ rdtice ! Mass flux by sublimation281 & * a_i_1d(ji) * r1_Dt_ice 282 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice ! Mass flux by sublimation 283 283 284 284 ! new snow thickness … … 337 337 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 338 338 339 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_ rdtice ! Heat flux to the ocean [W.m-2], <0339 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 340 340 ! ice enthalpy zEi is "sent" to the ocean 341 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_ rdtice ! Salt flux341 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 342 342 ! using s_i_1d and not sz_i_1d(jk) is ok 343 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_ rdtice ! Mass flux343 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 344 344 345 345 ELSE !-- Surface melting … … 363 363 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 364 364 365 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_ rdtice ! Salt flux >0365 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 366 366 ! using s_i_1d and not sz_i_1d(jk) is ok) 367 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_ rdtice ! Heat flux [W.m-2], < 0368 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_ rdtice ! Heat flux used in this process [W.m-2], > 0367 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux [W.m-2], < 0 368 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat flux used in this process [W.m-2], > 0 369 369 ! 370 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_ rdtice ! Mass flux370 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 371 371 372 372 END IF … … 378 378 dh_i_sub(ji) = dh_i_sub(ji) + zdum 379 379 380 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_ rdtice ! Salt flux >0380 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 381 381 ! clem: flux is sent to the ocean for simplicity 382 382 ! but salt should remain in the ice except 383 383 ! if all ice is melted. => must be corrected 384 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_ rdtice ! Heat flux [W.m-2], < 0385 386 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_ rdtice ! Mass flux > 0384 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0 385 386 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_Dt_ice ! Mass flux > 0 387 387 388 388 ! update remaining mass flux … … 409 409 ! remaining "potential" evap is sent to ocean 410 410 DO ji = 1, npti 411 wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_ rdtice ! <=0 (net evap for the ocean in kg.m-2.s-1)411 wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_Dt_ice ! <=0 (net evap for the ocean in kg.m-2.s-1) 412 412 END DO 413 413 … … 437 437 !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 438 438 !--- zswi2 if dh/dt > 3.6e-7 439 zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_ rdtice , epsi10 ) )439 zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_Dt_ice , epsi10 ) ) 440 440 zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 441 441 zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) … … 463 463 zfmdt = - rhoic * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0) 464 464 465 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_ rdtice ! Heat flux to the ocean [W.m-2], >0466 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_ rdtice ! Heat flux used in this process [W.m-2], <0467 468 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_ rdtice ! Salt flux, <0469 470 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * r1_ rdtice ! Mass flux, <0465 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux to the ocean [W.m-2], >0 466 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat flux used in this process [W.m-2], <0 467 468 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice ! Salt flux, <0 469 470 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice ! Mass flux, <0 471 471 472 472 ! update heat content (J.m-2) and layer thickness … … 499 499 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 500 500 501 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_ rdtice ! Heat flux to the ocean [W.m-2], <0501 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 502 502 ! ice enthalpy zEi is "sent" to the ocean 503 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_ rdtice ! Salt flux503 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 504 504 ! using s_i_1d and not sz_i_1d(jk) is ok 505 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_ rdtice ! Mass flux505 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 506 506 507 507 ! update heat content (J.m-2) and layer thickness … … 529 529 zQm = zfmdt * zEw ! Heat exchanged with ocean 530 530 531 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_ rdtice ! Heat flux to the ocean [W.m-2], <0532 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_ rdtice ! Heat used in this process [W.m-2], >0533 534 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_ rdtice ! Salt flux531 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 532 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat used in this process [W.m-2], >0 533 534 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 535 535 ! using s_i_1d and not sz_i_1d(jk) is ok 536 536 537 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_ rdtice ! Mass flux537 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 538 538 539 539 ! update heat content (J.m-2) and layer thickness … … 565 565 566 566 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 567 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_ rdtice ! Heat used to melt snow, W.m-2 (>0)568 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_ rdtice ! Mass flux567 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_Dt_ice ! Heat used to melt snow, W.m-2 (>0) 568 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! Mass flux 569 569 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 570 570 ! 571 571 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 572 hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_ rdtice572 hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 573 573 574 574 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) … … 580 580 ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 581 581 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 582 z1_rho = 1._wp / ( rhosn+r au0-rhoic )582 z1_rho = 1._wp / ( rhosn+rho0-rhoic ) 583 583 DO ji = 1, npti 584 584 ! 585 dh_snowice(ji) = MAX( 0._wp , ( rhosn * h_s_1d(ji) + (rhoic-r au0) * h_i_1d(ji) ) * z1_rho )585 dh_snowice(ji) = MAX( 0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rho0) * h_i_1d(ji) ) * z1_rho ) 586 586 587 587 h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji) … … 593 593 zQm = zfmdt * zEw 594 594 595 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_ rdtice ! Heat flux596 597 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_ rdtice ! Salt flux595 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux 596 597 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice ! Salt flux 598 598 599 599 ! Case constant salinity in time: virtual salt flux to keep salinity constant 600 600 IF( nn_icesal /= 2 ) THEN 601 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_ rdtice & ! put back sss_m into the ocean602 & - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_ rdtice ! and get rn_icesal from the ocean601 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice & ! put back sss_m into the ocean 602 & - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_Dt_ice ! and get rn_icesal from the ocean 603 603 ENDIF 604 604 605 605 ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 606 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_ rdtice607 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_ rdtice606 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_Dt_ice 607 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_Dt_ice 608 608 609 609 ! update heat content (J.m-2) and layer thickness -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_do.F90
r9604 r9923 140 140 ! Physical constants 141 141 zhicrit = 0.04 ! frazil ice thickness 142 ztwogp = 2. * r au0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav142 ztwogp = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoic ) ) ! reduced grav 143 143 zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag) 144 144 zgamafr = 0.03 … … 289 289 290 290 ! Contribution to heat flux to the ocean [W.m-2], >0 291 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_ rdtice291 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_Dt_ice 292 292 ! Total heat flux used in this process [W.m-2] 293 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_ rdtice293 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_Dt_ice 294 294 ! mass flux 295 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_ rdtice295 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_Dt_ice 296 296 ! salt flux 297 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_ rdtice297 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_Dt_ice 298 298 END DO 299 299 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_ent.F90
r9604 r9923 129 129 ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 130 130 DO ji = 1, npti 131 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_ rdtice * &131 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_Dt_ice * & 132 132 & ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 133 133 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_pnd.F90
r9750 r9923 168 168 ! melt pond mass flux (<0) 169 169 IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN 170 zfac = zfr_mlt * zdv_mlt * rhofw * r1_ rdtice170 zfac = zfr_mlt * zdv_mlt * rhofw * r1_Dt_ice 171 171 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 172 172 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_sal.F90
r9750 r9923 98 98 99 99 ! Salt flux 100 sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_ rdtice100 sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_Dt_ice 101 101 ENDIF 102 102 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_zdf_bl99.F90
r9656 r9923 770 770 771 771 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 772 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_ rdtice )*a_i_1d(ji)772 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice )*a_i_1d(ji) 773 773 ELSE ! case T_su = 0degC 774 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_ rdtice )*a_i_1d(ji)774 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice )*a_i_1d(ji) 775 775 ENDIF 776 776 777 777 ELSEIF( k_jules == np_jules_ACTIVE ) THEN 778 778 779 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_ rdtice ) * a_i_1d(ji)779 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice ) * a_i_1d(ji) 780 780 781 781 ENDIF … … 785 785 ! 786 786 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 787 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_ rdtice * a_i_1d(ji)787 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_Dt_ice * a_i_1d(ji) 788 788 ! 789 789 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceupdate.F90
r9784 r9923 174 174 snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) 175 175 ! ! time evolution of snow+ice mass 176 snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_ rdtice176 snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice 177 177 178 178 END DO … … 336 336 ENDIF 337 337 338 zrhoco = r au0 * rn_cio338 zrhoco = rho0 * rn_cio 339 339 ! 340 340 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icevar.F90
r9725 r9923 477 477 DO ji = 1 , jpi 478 478 ! update exchanges with ocean 479 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_ rdtice ! W.m-2 <0479 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 480 480 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 481 481 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) … … 488 488 DO ji = 1 , jpi 489 489 ! update exchanges with ocean 490 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_ rdtice ! W.m-2 <0490 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 491 491 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj) 492 492 t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) … … 498 498 DO ji = 1 , jpi 499 499 ! update exchanges with ocean 500 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoic * r1_ rdtice501 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_ rdtice502 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_ rdtice500 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoic * r1_Dt_ice 501 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_Dt_ice 502 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_Dt_ice 503 503 ! 504 504 !----------------------------------------------------------------- … … 669 669 ! In case snow load is in excess that would lead to transformation from snow to ice 670 670 ! Then, transfer the snow excess into the ice (different from icethd_dh) 671 zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - r au0 ) * zh_i(ji,jl) ) * r1_rau0 )671 zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rho0 ) * zh_i(ji,jl) ) * r1_rho0 ) 672 672 ! recompute h_i, h_s avoiding out of bounds values 673 673 zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icewri.F90
r9604 r9923 85 85 ! Standard outputs 86 86 !----------------- 87 zrho1 = ( r au0 - rhoic ) * r1_rau0; zrho2 = rhosn * r1_rau087 zrho1 = ( rho0 - rhoic ) * r1_rho0 ; zrho2 = rhosn * r1_rho0 88 88 ! masks 89 89 IF( iom_use('icemask' ) ) CALL iom_put( "icemask" , zmsk00 ) ! ice mask 0% … … 250 250 CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up") 251 251 252 CALL histdef( kid, "sithic", "Ice thickness" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )253 CALL histdef( kid, "siconc", "Ice concentration" , "%" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )254 CALL histdef( kid, "sitemp", "Ice temperature" , "C" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )255 CALL histdef( kid, "sivelu", "i-Ice speed " , "m/s" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )256 CALL histdef( kid, "sivelv", "j-Ice speed " , "m/s" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )257 CALL histdef( kid, "sistru", "i-Wind stress over ice" , "Pa" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )258 CALL histdef( kid, "sistrv", "j-Wind stress over ice" , "Pa" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )259 CALL histdef( kid, "sisflx", "Solar flx over ocean" , "W/m2" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )260 CALL histdef( kid, "sinflx", "NonSolar flx over ocean", "W/m2" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )261 CALL histdef( kid, "snwpre", "Snow precipitation" , "kg/m2/s", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )262 CALL histdef( kid, "sisali", "Ice salinity" , "PSU" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )263 CALL histdef( kid, "sivolu", "Ice volume" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )264 CALL histdef( kid, "sidive", "Ice divergence" , "10-8s-1", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )265 CALL histdef( kid, "si_amp", "Melt pond fraction" , "%" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )266 CALL histdef( kid, "si_vmp", "Melt pond volume" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", r dt, rdt )267 ! 268 CALL histdef( kid, "sithicat", "Ice thickness" , "m" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )269 CALL histdef( kid, "siconcat", "Ice concentration" , "%" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )270 CALL histdef( kid, "sisalcat", "Ice salinity" , "" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )271 CALL histdef( kid, "snthicat", "Snw thickness" , "m" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )252 CALL histdef( kid, "sithic", "Ice thickness" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 253 CALL histdef( kid, "siconc", "Ice concentration" , "%" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 254 CALL histdef( kid, "sitemp", "Ice temperature" , "C" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 255 CALL histdef( kid, "sivelu", "i-Ice speed " , "m/s" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 256 CALL histdef( kid, "sivelv", "j-Ice speed " , "m/s" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 257 CALL histdef( kid, "sistru", "i-Wind stress over ice" , "Pa" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 258 CALL histdef( kid, "sistrv", "j-Wind stress over ice" , "Pa" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 259 CALL histdef( kid, "sisflx", "Solar flx over ocean" , "W/m2" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 260 CALL histdef( kid, "sinflx", "NonSolar flx over ocean", "W/m2" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 261 CALL histdef( kid, "snwpre", "Snow precipitation" , "kg/m2/s", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 262 CALL histdef( kid, "sisali", "Ice salinity" , "PSU" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 263 CALL histdef( kid, "sivolu", "Ice volume" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 264 CALL histdef( kid, "sidive", "Ice divergence" , "10-8s-1", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 265 CALL histdef( kid, "si_amp", "Melt pond fraction" , "%" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 266 CALL histdef( kid, "si_vmp", "Melt pond volume" , "m" , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 267 ! 268 CALL histdef( kid, "sithicat", "Ice thickness" , "m" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rn_Dt, rn_Dt ) 269 CALL histdef( kid, "siconcat", "Ice concentration" , "%" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rn_Dt, rn_Dt ) 270 CALL histdef( kid, "sisalcat", "Ice salinity" , "" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rn_Dt, rn_Dt ) 271 CALL histdef( kid, "snthicat", "Snw thickness" , "m" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rn_Dt, rn_Dt ) 272 272 273 273 CALL histend( kid, snc4set ) ! end of the file definition -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/NST/agrif_oce_update.F90
r9863 r9923 381 381 END DO 382 382 383 IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 383 IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part 384 384 385 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 385 ! Add asselin part386 386 DO jn = n1, n2-1 387 387 DO jk = 1, jpk … … 389 389 DO ji = i1, i2 390 390 IF( tabres_child(ji,jj,jk,jn) /= 0. ) THEN 391 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 392 & + atfp * ( tabres_child(ji,jj,jk,jn) & 393 & - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 391 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 392 & + rn_atfp * ( tabres_child(ji,jj,jk,jn) - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 394 393 ENDIF 395 394 END DO … … 458 457 ztnu = tabres(ji,jj,jk,jn) 459 458 ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 460 tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 461 & * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 459 tsb(ji,jj,jk,jn) = ( ztb + rn_rn_atfp * ( ztnu - ztno) ) / e3t_b(ji,jj,jk) * tmask(ji,jj,jk) 462 460 ENDIF 463 461 END DO … … 568 566 END DO 569 567 570 DO jk =1,jpk571 DO jj =j1,j2572 DO ji =i1,i2568 DO jk = 1, jpk 569 DO jj = j1, j2 570 DO ji = i1, i2 573 571 IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler) ) THEN ! Add asselin part 574 572 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 575 ub(ji,jj,jk) = ub(ji,jj,jk) & 576 & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 573 ub(ji,jj,jk) = ub(ji,jj,jk) + rn_atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 577 574 ENDIF 578 !579 575 un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 580 576 END DO … … 615 611 zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk) 616 612 zunu = tabres(ji,jj,jk,1) 617 ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) & 618 & * umask(ji,jj,jk) / e3u_b(ji,jj,jk) 613 ub(ji,jj,jk) = ( zub + rn_atfp * ( zunu - zuno) ) / e3u_b(ji,jj,jk) * umask(ji,jj,jk) 619 614 ENDIF 620 615 ! 621 un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)616 un(ji,jj,jk) = tabres(ji,jj,jk,1) / e3u_n(ji,jj,jk) * umask(ji,jj,jk) 622 617 END DO 623 618 END DO … … 761 756 IF( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part 762 757 !!gm IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 763 vb(ji,jj,jk) = vb(ji,jj,jk) & 764 & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 758 vb(ji,jj,jk) = vb(ji,jj,jk) + rn_atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 765 759 ENDIF 766 760 ! … … 807 801 zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk) 808 802 zvnu = tabres(ji,jj,jk,1) 809 vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) & 810 & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk) 803 vb(ji,jj,jk) = ( zvb + rn_atfp * ( zvnu - zvno) ) / e3v_b(ji,jj,jk) * vmask(ji,jj,jk) 811 804 ENDIF 812 805 ! 813 vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)806 vn(ji,jj,jk) = tabres(ji,jj,jk,1) / e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 814 807 END DO 815 808 END DO … … 913 906 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 914 907 zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj) 915 ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)908 ub_b(ji,jj) = ub_b(ji,jj) + rn_atfp * zcorr * umask(ji,jj,1) 916 909 END IF 917 910 ENDIF … … 981 974 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 982 975 zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj) 983 vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)976 vb_b(ji,jj) = vb_b(ji,jj) + rn_atfp * zcorr * vmask(ji,jj,1) 984 977 END IF 985 978 ENDIF … … 1032 1025 DO jj=j1,j2 1033 1026 DO ji=i1,i2 1034 sshb(ji,jj) = sshb(ji,jj) & 1035 & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 1027 sshb(ji,jj) = sshb(ji,jj) + rn_atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 1036 1028 END DO 1037 1029 END DO … … 1126 1118 IF ( western_side ) THEN 1127 1119 DO jj=j1,j2 1128 zcor = r dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))1120 zcor = rn_Dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 1129 1121 sshn(i1 ,jj) = sshn(i1 ,jj) + zcor 1130 IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) sshb(i1 ,jj) = sshb(i1 ,jj) + atfp * zcor1131 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1 ,jj) = sshb(i1 ,jj) + atfp * zcor1122 IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) sshb(i1 ,jj) = sshb(i1 ,jj) + rn_atfp * zcor 1123 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1 ,jj) = sshb(i1 ,jj) + rn_atfp * zcor 1132 1124 END DO 1133 1125 ENDIF 1134 1126 IF ( eastern_side ) THEN 1135 1127 DO jj=j1,j2 1136 zcor = - r dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))1128 zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 1137 1129 sshn(i2+1,jj) = sshn(i2+1,jj) + zcor 1138 IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor1139 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor1130 IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) sshb(i2+1,jj) = sshb(i2+1,jj) + rn_atfp * zcor 1131 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + rn_atfp * zcor 1140 1132 END DO 1141 1133 ENDIF … … 1218 1210 IF (southern_side) THEN 1219 1211 DO ji=i1,i2 1220 zcor = r dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1))1212 zcor = rn_Dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1)) 1221 1213 sshn(ji,j1 ) = sshn(ji,j1 ) + zcor 1222 IF ( .NOT.( lk_agrif_fstep .AND. l_euler ) ) sshb(ji,j1 ) = sshb(ji,j1) + atfp * zcor1223 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1 ) = sshb(ji,j1) + atfp * zcor1214 IF ( .NOT.( lk_agrif_fstep .AND. l_euler ) ) sshb(ji,j1 ) = sshb(ji,j1) + rn_atfp * zcor 1215 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1 ) = sshb(ji,j1) + rn_atfp * zcor 1224 1216 END DO 1225 1217 ENDIF 1226 1218 IF (northern_side) THEN 1227 1219 DO ji=i1,i2 1228 zcor = - r dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2))1220 zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2)) 1229 1221 sshn(ji,j2+1) = sshn(ji,j2+1) + zcor 1230 IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor1231 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor1222 IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) sshb(ji,j2+1) = sshb(ji,j2+1) + rn_atfp * zcor 1223 !!gm IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + rn_atfp * zcor 1232 1224 END DO 1233 1225 ENDIF … … 1381 1373 DO jj=j1,j2 1382 1374 DO ji=i1,i2 1383 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) & 1384 & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) ) 1375 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) + rn_atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) ) 1385 1376 END DO 1386 1377 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/NST/agrif_top_update.F90
r9863 r9923 140 140 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 141 141 trb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 142 & + atfp * ( tabres_child(ji,jj,jk,jn) & 143 & - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 142 & + rn_atfp * ( tabres_child(ji,jj,jk,jn) - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 144 143 ENDIF 145 144 END DO … … 209 208 ztnu = tabres(ji,jj,jk,jn) 210 209 ztno = trn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 211 trb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 212 & * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 210 trb(ji,jj,jk,jn) = ( ztb + rn_atfp * ( ztnu - ztno) ) / e3t_b(ji,jj,jk) * tmask(ji,jj,jk) 213 211 ENDIF 214 212 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/NST/agrif_user.F90
r9788 r9923 217 217 218 218 ! Check time steps 219 IF( NINT(Agrif_Rhot()) * NINT(r dt) .NE. Agrif_Parent(rdt) ) THEN220 WRITE(cl_check1,*) NINT(Agrif_Parent(r dt))221 WRITE(cl_check2,*) NINT(r dt)222 WRITE(cl_check3,*) NINT(Agrif_Parent(r dt)/Agrif_Rhot())219 IF( NINT(Agrif_Rhot()) * NINT(rn_Dt) /= Agrif_Parent(rn_Dt) ) THEN 220 WRITE(cl_check1,*) NINT(Agrif_Parent(rn_Dt)) 221 WRITE(cl_check2,*) NINT(rn_Dt) 222 WRITE(cl_check3,*) NINT(Agrif_Parent(rn_Dt)/Agrif_Rhot()) 223 223 CALL ctl_stop( 'Incompatible time step between ocean grids', & 224 224 & 'parent grid value : '//cl_check1 , & … … 229 229 ! Check run length 230 230 IF( Agrif_IRhot() * (Agrif_Parent(nitend)- & 231 Agrif_Parent(nit000)+1) .NE.(nitend-nit000+1) ) THEN231 Agrif_Parent(nit000)+1) /= (nitend-nit000+1) ) THEN 232 232 WRITE(cl_check1,*) (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 233 233 WRITE(cl_check2,*) Agrif_Parent(nitend) *Agrif_IRhot() … … 601 601 IF( check_namelist ) THEN 602 602 ! Check time steps 603 IF( NINT(Agrif_Rhot()) * NINT(r dt) .NE. Agrif_Parent(rdt) ) THEN604 WRITE(cl_check1,*) Agrif_Parent(r dt)605 WRITE(cl_check2,*) r dt606 WRITE(cl_check3,*) r dt*Agrif_Rhot()603 IF( NINT(Agrif_Rhot()) * NINT(rn_Dt) .NE. Agrif_Parent(rn_Dt) ) THEN 604 WRITE(cl_check1,*) Agrif_Parent(rn_Dt) 605 WRITE(cl_check2,*) rn_Dt 606 WRITE(cl_check3,*) rn_Dt*Agrif_Rhot() 607 607 CALL ctl_stop( 'incompatible time step between grids', & 608 608 & 'parent grid value : '//cl_check1 , & -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/ASM/asminc.F90
r9863 r9923 536 536 ! 537 537 it = kt - nit000 + 1 538 zincwgt = wgtiau(it) / r dt ! IAU weight for the current time step538 zincwgt = wgtiau(it) / rn_Dt ! IAU weight for the current time step 539 539 ! 540 540 IF(lwp) THEN … … 651 651 ! 652 652 it = kt - nit000 + 1 653 zincwgt = wgtiau(it) / r dt ! IAU weight for the current time step653 zincwgt = wgtiau(it) / rn_Dt ! IAU weight for the current time step 654 654 ! 655 655 IF(lwp) THEN … … 721 721 ! 722 722 it = kt - nit000 + 1 723 zincwgt = wgtiau(it) / r dt ! IAU weight for the current time step723 zincwgt = wgtiau(it) / rn_Dt ! IAU weight for the current time step 724 724 ! 725 725 IF(lwp) THEN … … 841 841 it = kt - nit000 + 1 842 842 zincwgt = wgtiau(it) ! IAU weight for the current time step 843 ! note this is not a tendency so should not be divided by r dt (as with the tracer and other increments)843 ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments) 844 844 ! 845 845 IF(lwp) THEN … … 876 876 #if defined key_cice && defined key_asminc 877 877 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 878 ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / r dt878 ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rn_Dt 879 879 #endif 880 880 ! … … 926 926 #if defined key_cice && defined key_asminc 927 927 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 928 ndaice_da(:,:) = seaice_bkginc(:,:) / r dt928 ndaice_da(:,:) = seaice_bkginc(:,:) / rn_Dt 929 929 #endif 930 930 IF ( .NOT. PRESENT(kindic) ) THEN … … 959 959 ! ! fwf : ice formation and melting 960 960 ! 961 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) ) *rdt961 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) ) * rn_Dt 962 962 ! 963 963 ! ! change salinity down to mixed layer depth … … 1008 1008 ! !! ! E-P (kg m-2 s-2) 1009 1009 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) 1010 ! END DO !ji1011 ! END DO !jj!1010 ! END DO !ji 1011 ! END DO !jj! 1012 1012 ! 1013 1013 ! ENDIF !ln_seaicebal -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/BDY/bdyice.F90
r9657 r9923 124 124 125 125 ! Then, a) transfer the snow excess into the ice (different from icethd_dh) 126 zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - r au0 ) * h_i(ji,jj,jl) ) * r1_rau0 )126 zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 ) 127 127 ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment) 128 128 !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhosn / rhoic ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/BDY/bdylib.F90
r9598 r9923 4 4 !! Unstructured Open Boundary Cond. : Library module of generic boundary algorithms. 5 5 !!====================================================================== 6 !! History : 3.6 ! 2013 (D. Storkey) original code7 !! 4.0 ! 2014 (T. Lovato) Generalize OBC structure6 !! History : 3.6 ! 2013 (D. Storkey) original code 7 !! 4.0 ! 2014 (T. Lovato) Generalize OBC structure 8 8 !!---------------------------------------------------------------------- 9 9 10 !!---------------------------------------------------------------------- 10 !! bdy_orlanski_2d 11 !! bdy_orlanski_3d 11 !! bdy_frs : Apply the Flow Relaxation Scheme (tracers) 12 !! bdy_spe : Apply a specified value (tracers) 13 !! bdy_orl : Apply Orlanski radiation (tracers) 14 !! bdy_orlanski_2d: 2D - - - 15 !! bdy_orlanski_3d: 3D - - - 16 !! bdy_nmn : Duplicate the value at open boundaries (zero gradient) 12 17 !!---------------------------------------------------------------------- 13 18 USE oce ! ocean dynamics and tracers … … 22 27 PRIVATE 23 28 24 PUBLIC bdy_frs, bdy_spe, bdy_nmn, bdy_orl 25 PUBLIC bdy_orlanski_2d 26 PUBLIC bdy_orlanski_3d 29 PUBLIC bdy_frs, bdy_spe, bdy_nmn 30 PUBLIC bdy_orl, bdy_orlanski_2d, bdy_orlanski_3d 27 31 28 32 !!---------------------------------------------------------------------- … … 230 234 ! Note no rdt factor in expression for zdt because it cancels in the expressions for 231 235 ! zrx and zry. 232 zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1)236 zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1) 233 237 zdx = ( ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) ) / zex2 ) * zmask_x 234 238 zdy_1 = ( ( phib(iibm1 ,ijbm1 ) - phib(iibm1jm1,ijbm1jm1) ) / zey1 ) * zmask_y1 … … 247 251 zout = sign( 1., zrx ) 248 252 zout = 0.5*( zout + abs(zout) ) 249 zwgt = 2.*rdt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )253 zwgt = rDt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) ) 250 254 ! only apply radiation on outflow points 251 255 if( ll_npo ) then !! NPO version !! … … 385 389 ! Centred derivative is calculated as average of "left" and "right" derivatives for 386 390 ! this reason. 387 zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk)391 zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk) 388 392 zdx = ( ( phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk) ) / zex2 ) * zmask_x 389 393 zdy_1 = ( ( phib(iibm1 ,ijbm1 ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) / zey1 ) * zmask_y1 … … 403 407 zout = sign( 1., zrx ) 404 408 zout = 0.5*( zout + abs(zout) ) 405 zwgt = 2.*rdt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) )409 zwgt = rDt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) ) 406 410 ! only apply radiation on outflow points 407 411 if( ll_npo ) then !! NPO version !! … … 426 430 ! 427 431 END SUBROUTINE bdy_orlanski_3d 432 428 433 429 434 SUBROUTINE bdy_nmn( idx, igrd, phia ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/BDY/bdytides.F90
r9598 r9923 295 295 !!---------------------------------------------------------------------- 296 296 ! 297 ilen0(1) = SIZE( td%ssh(:,1,1))298 ilen0(2) = SIZE( td%u(:,1,1))299 ilen0(3) = SIZE( td%v(:,1,1))297 ilen0(1) = SIZE( td%ssh(:,1,1) ) 298 ilen0(2) = SIZE( td%u (:,1,1) ) 299 ilen0(3) = SIZE( td%v (:,1,1) ) 300 300 301 301 zflag=1 302 302 IF ( PRESENT(jit) ) THEN 303 IF ( jit /= 1 ) zflag=0303 IF ( jit /= 1 ) zflag=0 304 304 ENDIF 305 305 306 IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN306 IF ( ( nsec_day == NINT( 0.5_wp * rn_Dt ) .OR. kt == nit000 ) .AND. zflag==1 ) THEN 307 307 ! 308 kt_tide = kt - (nsec_day - 0.5_wp * r dt)/rdt308 kt_tide = kt - (nsec_day - 0.5_wp * rn_Dt) / rn_Dt 309 309 ! 310 310 IF(lwp) THEN 311 311 WRITE(numout,*) 312 WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=', kt312 WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=', kt 313 313 WRITE(numout,*) '~~~~~~~~~~~~~~ ' 314 314 ENDIF … … 325 325 326 326 IF( PRESENT(jit) ) THEN 327 z_arg = ((kt-kt_tide) * r dt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )327 z_arg = ((kt-kt_tide) * rn_Dt + (jit+0.5_wp*(time_add-1)) * rn_Dt / REAL(nn_e,wp) ) 328 328 ELSE 329 z_arg = ((kt-kt_tide)+time_add) * r dt329 z_arg = ((kt-kt_tide)+time_add) * rn_Dt 330 330 ENDIF 331 331 332 332 ! Linear ramp on tidal component at open boundaries 333 333 zramp = 1._wp 334 IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*r dt)/(rdttideramp*rday),0._wp),1._wp)334 IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rn_Dt)/(rn_ramp*rday),0._wp),1._wp) 335 335 336 336 DO itide = 1, nb_harmo … … 392 392 ! Absolute time from model initialization: 393 393 IF( PRESENT(kit) ) THEN 394 z_arg = ( kt + (kit+time_add-1) / REAL(nn_ baro,wp) ) * rdt394 z_arg = ( kt + (kit+time_add-1) / REAL(nn_e,wp) ) * rn_Dt 395 395 ELSE 396 z_arg = ( kt + time_add ) * r dt396 z_arg = ( kt + time_add ) * rn_Dt 397 397 ENDIF 398 398 399 399 ! Linear ramp on tidal component at open boundaries 400 400 zramp = 1. 401 IF ( ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)401 IF ( ln_tide_ramp ) zramp = MIN( MAX( 0. , (z_arg - nit000*rn_Dt)/(rn_ramp*rday) ) , 1. ) 402 402 403 403 DO ib_bdy = 1,nb_bdy … … 414 414 ! We refresh nodal factors every day below 415 415 ! This should be done somewhere else 416 IF ( ( nsec_day == NINT(0.5_wp * r dt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN417 ! 418 kt_tide = kt - (nsec_day - 0.5_wp * r dt)/rdt416 IF ( ( nsec_day == NINT(0.5_wp * rn_Dt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN 417 ! 418 kt_tide = kt - (nsec_day - 0.5_wp * rn_Dt) / rn_Dt 419 419 ! 420 420 IF(lwp) THEN … … 428 428 ! 429 429 ENDIF 430 zoff = -kt_tide * r dt! time offset relative to nodal factor computation time430 zoff = -kt_tide * rn_Dt ! time offset relative to nodal factor computation time 431 431 ! 432 432 ! If time splitting, initialize arrays from slow varying open boundary data: -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/BDY/bdyvol.F90
r9598 r9923 84 84 ! ----------------------------------------------------------------------- 85 85 !!gm replace these lines : 86 z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / r au086 z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rho0 87 87 IF( lk_mpp ) CALL mpp_sum( z_cflxemp ) ! sum over the global domain 88 88 !!gm by : 89 !!gm z_cflxemp = glob_sum( ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / r au089 !!gm z_cflxemp = glob_sum( ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rho0 90 90 !!gm 91 91 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/dia25h.F90
r9598 r9923 139 139 ! ----------------- 140 140 ! Define frequency of summing to create 25 h mean 141 IF( MOD( 3600 ,INT(rdt) ) == 0 ) THEN142 i_steps = 3600 /INT(rdt)141 IF( MOD( 3600 , INT(rn_Dt) ) == 0 ) THEN 142 i_steps = 3600 / INT( rn_Dt ) 143 143 ELSE 144 CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,r dt) = 0 otherwise no hourly values are possible')144 CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rn_Dt) = 0 otherwise no hourly values are possible') 145 145 ENDIF 146 146 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diaar5.F90
r9598 r9923 161 161 162 162 ! ! ocean bottom pressure 163 zztmp = r au0 * grav * 1.e-4_wp ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa163 zztmp = rho0 * grav * 1.e-4_wp ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 164 164 zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) ) 165 165 CALL iom_put( 'botpres', zbotpres ) … … 198 198 END IF 199 199 ! 200 zmass = r au0 * ( zarho + zvol ) ! total mass of liquid seawater200 zmass = rho0 * ( zarho + zvol ) ! total mass of liquid seawater 201 201 ztemp = ztemp / zvol ! potential temperature in liquid seawater 202 202 zsal = zsal / zvol ! Salinity of liquid seawater … … 239 239 DO ji = 1, jpi 240 240 DO jj = 1, jpj 241 zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * r au0 * e3w_n(ji, jj, jk)241 zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rho0 * e3w_n(ji, jj, jk) 242 242 END DO 243 243 END DO … … 287 287 CALL lbc_lnk( z2d, 'U', -1. ) 288 288 IF( cptr == 'adv' ) THEN 289 IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , r au0_rcp * z2d ) ! advective heat transport in i-direction290 IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , r au0 * z2d ) ! advective salt transport in i-direction289 IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , rho0_rcp * z2d ) ! advective heat transport in i-direction 290 IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , rho0 * z2d ) ! advective salt transport in i-direction 291 291 ENDIF 292 292 IF( cptr == 'ldf' ) THEN 293 IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , r au0_rcp * z2d ) ! diffusive heat transport in i-direction294 IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , r au0 * z2d ) ! diffusive salt transport in i-direction293 IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , rho0_rcp * z2d ) ! diffusive heat transport in i-direction 294 IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , rho0 * z2d ) ! diffusive salt transport in i-direction 295 295 ENDIF 296 296 ! … … 305 305 CALL lbc_lnk( z2d, 'V', -1. ) 306 306 IF( cptr == 'adv' ) THEN 307 IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , r au0_rcp * z2d ) ! advective heat transport in j-direction308 IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , r au0 * z2d ) ! advective salt transport in j-direction307 IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , rho0_rcp * z2d ) ! advective heat transport in j-direction 308 IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , rho0 * z2d ) ! advective salt transport in j-direction 309 309 ENDIF 310 310 IF( cptr == 'ldf' ) THEN 311 IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , r au0_rcp * z2d ) ! diffusive heat transport in j-direction312 IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , r au0 * z2d ) ! diffusive salt transport in j-direction311 IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , rho0_rcp * z2d ) ! diffusive heat transport in j-direction 312 IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , rho0 * z2d ) ! diffusive salt transport in j-direction 313 313 ENDIF 314 314 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diacfl.F90
r9863 r9923 66 66 DO jj = 1, jpj 67 67 DO ji = 1, fs_jpim1 ! vector opt. 68 zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * r 2dt / e1u (ji,jj) ! for i-direction69 zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * r 2dt / e2v (ji,jj) ! for j-direction70 zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * r 2dt / e3w_n(ji,jj,jk) ! for k-direction68 zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * rDt / e1u (ji,jj) ! for i-direction 69 zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * rDt / e2v (ji,jj) ! for j-direction 70 zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * rDt / e3w_n(ji,jj,jk) ! for k-direction 71 71 END DO 72 72 END DO … … 115 115 WRITE(numcfl,*) '******************************************' 116 116 WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3) 117 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r 2dt/rCu_max117 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCu_max 118 118 WRITE(numcfl,*) '******************************************' 119 119 WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3) 120 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r 2dt/rCv_max120 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCv_max 121 121 WRITE(numcfl,*) '******************************************' 122 122 WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3) 123 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r 2dt/rCw_max123 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCw_max 124 124 CLOSE( numcfl ) 125 125 ! … … 128 128 WRITE(numout,*) 'dia_cfl : Maximum Courant number information for the run ' 129 129 WRITE(numout,*) '~~~~~~~' 130 WRITE(numout,*) ' Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', r 2dt/rCu_max131 WRITE(numout,*) ' Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', r 2dt/rCv_max132 WRITE(numout,*) ' Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', r 2dt/rCw_max130 WRITE(numout,*) ' Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', rDt/rCu_max 131 WRITE(numout,*) ' Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', rDt/rCv_max 132 WRITE(numout,*) ' Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', rDt/rCw_max 133 133 ENDIF 134 134 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diadct.F90
r9598 r9923 679 679 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 680 680 zrhop = interp(k%I,k%J,jk,'V',rhop) 681 zrhoi = interp(k%I,k%J,jk,'V',rhd*r au0+rau0)681 zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0) 682 682 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 683 683 CASE(2,3) … … 685 685 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 686 686 zrhop = interp(k%I,k%J,jk,'U',rhop) 687 zrhoi = interp(k%I,k%J,jk,'U',rhd*r au0+rau0)687 zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0) 688 688 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 689 689 END SELECT … … 851 851 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 852 852 zrhop = interp(k%I,k%J,jk,'V',rhop) 853 zrhoi = interp(k%I,k%J,jk,'V',rhd*r au0+rau0)853 zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0) 854 854 855 855 CASE(2,3) … … 857 857 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 858 858 zrhop = interp(k%I,k%J,jk,'U',rhop) 859 zrhoi = interp(k%I,k%J,jk,'U',rhd*r au0+rau0)859 zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0) 860 860 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 861 861 END SELECT -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diaharm.F90
r9598 r9923 181 181 IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN 182 182 ! 183 ztime = ( kt-nit000+1) * rdt183 ztime = ( kt - nit000+1 ) * rn_Dt 184 184 ! 185 185 nhc = 0 … … 231 231 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 232 232 233 ztime_ini = nit000_han*r dt ! Initial time in seconds at the beginning of analysis234 ztime_end = nitend_han*r dt ! Final time in seconds at the end of analysis233 ztime_ini = nit000_han*rn_Dt ! Initial time in seconds at the beginning of analysis 234 ztime_end = nitend_han*rn_Dt ! Final time in seconds at the end of analysis 235 235 nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis 236 236 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diahsb.F90
r9598 r9923 91 91 ! 1 - Trends due to forcing ! 92 92 ! ------------------------- ! 93 z_frc_trd_v = r1_r au0 * glob_sum( - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) ) ! volume fluxes93 z_frc_trd_v = r1_rho0 * glob_sum( - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) ) ! volume fluxes 94 94 z_frc_trd_t = glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) ) ! heat fluxes 95 95 z_frc_trd_s = glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) ) ! salt fluxes … … 100 100 IF( ln_isf ) z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) ) 101 101 ! ! Add penetrative solar radiation 102 IF( ln_traqsr ) z_frc_trd_t = z_frc_trd_t + r1_r au0_rcp * glob_sum( qsr (:,:) * surf(:,:) )102 IF( ln_traqsr ) z_frc_trd_t = z_frc_trd_t + r1_rho0_rcp * glob_sum( qsr (:,:) * surf(:,:) ) 103 103 ! ! Add geothermal heat flux 104 104 IF( ln_trabbc ) z_frc_trd_t = z_frc_trd_t + glob_sum( qgh_trd0(:,:) * surf(:,:) ) … … 120 120 ENDIF 121 121 122 frc_v = frc_v + z_frc_trd_v * r dt123 frc_t = frc_t + z_frc_trd_t * r dt124 frc_s = frc_s + z_frc_trd_s * r dt122 frc_v = frc_v + z_frc_trd_v * rn_Dt 123 frc_t = frc_t + z_frc_trd_t * rn_Dt 124 frc_s = frc_s + z_frc_trd_s * rn_Dt 125 125 ! ! Advection flux through fixed surface (z=0) 126 126 IF( ln_linssh ) THEN 127 frc_wn_t = frc_wn_t + z_wn_trd_t * r dt128 frc_wn_s = frc_wn_s + z_wn_trd_s * r dt127 frc_wn_t = frc_wn_t + z_wn_trd_t * rn_Dt 128 frc_wn_s = frc_wn_s + z_wn_trd_s * rn_Dt 129 129 ENDIF 130 130 … … 196 196 197 197 CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 ) ! vol - surface forcing (km3) 198 CALL iom_put( 'bgfrctem' , frc_t * r au0 * rcp * 1.e-20 )! hc - surface forcing (1.e20 J)199 CALL iom_put( 'bgfrchfx' , frc_t * r au0 * rcp / &! hc - surface forcing (W/m2)200 & ( surf_tot * kt * rdt ))198 CALL iom_put( 'bgfrctem' , frc_t * rho0_rcp * 1.e-20 ) ! hc - surface forcing (1.e20 J) 199 CALL iom_put( 'bgfrchfx' , frc_t * rho0_rcp / & ! hc - surface forcing (W/m2) 200 & ( surf_tot * kt * rn_Dt ) ) 201 201 CALL iom_put( 'bgfrcsal' , frc_s * 1.e-9 ) ! sc - surface forcing (psu*km3) 202 202 … … 204 204 CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot ) ! Temperature drift (C) 205 205 CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot ) ! Salinity drift (PSU) 206 CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * r au0 * rcp )! Heat content drift (1.e20 J)207 CALL iom_put( 'bgheatfx' , zdiff_hc * r au0 * rcp / &! Heat flux drift (W/m2)208 & ( surf_tot * kt * r dt ))206 CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rho0_rcp ) ! Heat content drift (1.e20 J) 207 CALL iom_put( 'bgheatfx' , zdiff_hc * rho0_rcp / & ! Heat flux drift (W/m2) 208 & ( surf_tot * kt * rn_Dt ) ) 209 209 CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9 ) ! Salt content drift (psu*km3) 210 210 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 ) ! volume ssh drift (km3) … … 224 224 CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot) ! Heat content drift (C) 225 225 CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot) ! Salt content drift (PSU) 226 CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * r au0 * rcp )! Heat content drift (1.e20 J)227 CALL iom_put( 'bgheatfx' , zdiff_hc1 * r au0 * rcp / &! Heat flux drift (W/m2)228 & ( surf_tot * kt * r dt ))226 CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rho0_rcp ) ! Heat content drift (1.e20 J) 227 CALL iom_put( 'bgheatfx' , zdiff_hc1 * rho0_rcp / & ! Heat flux drift (W/m2) 228 & ( surf_tot * kt * rn_Dt ) ) 229 229 CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9 ) ! Salt content drift (psu*km3) 230 230 CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 ) ! volume ssh drift (km3) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diahth.F90
r9598 r9923 89 89 REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth 90 90 REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth 91 REAL(wp) :: zthick_0 , zcoef ! temporaryscalars92 REAL(wp) :: zztmp, zzdep ! temporaryscalars inside do loop93 REAL(wp) :: zu, zv, zw, zut, zvt ! temporaryworkspace91 REAL(wp) :: zthick_0 ! local scalars 92 REAL(wp) :: zztmp, zzdep ! local scalars inside do loop 93 REAL(wp) :: zu, zv, zw, zut, zvt ! local workspace 94 94 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2 95 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2 … … 328 328 END DO 329 329 ! from temperature to heat contain 330 zcoef = rau0 * rcp 331 htc3(:,:) = zcoef * htc3(:,:) 330 htc3(:,:) = rho0_rcp * htc3(:,:) 332 331 CALL iom_put( "hc300", htc3 ) ! first 300m heat content 333 332 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/dianam.F90
r9598 r9923 71 71 ENDIF 72 72 73 IF( llfsec .OR. kfreq < 0 ) THEN ; inbsec = kfreq 74 ELSE ; inbsec = kfreq * NINT( r dt )! from time-step to seconds73 IF( llfsec .OR. kfreq < 0 ) THEN ; inbsec = kfreq ! output frequency already in seconds 74 ELSE ; inbsec = kfreq * NINT( rn_Dt ) ! from time-step to seconds 75 75 ENDIF 76 76 iddss = NINT( rday ) ! number of seconds in 1 day … … 116 116 ! date of the beginning and the end of the run 117 117 118 zdrun = r dt / rday * REAL( nitend - nit000, wp )! length of the run in days119 zjul = fjulday - r dt / rday118 zdrun = rn_Dt / rday * REAL( nitend - nit000, wp ) ! length of the run in days 119 zjul = fjulday - rn_Dt / rday 120 120 CALL ju2ymds( zjul , iyear1, imonth1, iday1, zsec1 ) ! year/month/day of the beginning of run 121 121 CALL ju2ymds( zjul + zdrun, iyear2, imonth2, iday2, zsec2 ) ! year/month/day of the end of run -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diaptr.F90
r9598 r9923 52 52 53 53 REAL(wp) :: rc_sv = 1.e-6_wp ! conversion from m3/s to Sverdrup 54 REAL(wp) :: rc_pwatt = 1.e-15_wp ! conversion from W to PW (further x r au0 x Cp)54 REAL(wp) :: rc_pwatt = 1.e-15_wp ! conversion from W to PW (further x rho0 x Cp) 55 55 REAL(wp) :: rc_ggram = 1.e-6_wp ! conversion from g to Pg 56 56 … … 424 424 IF( dia_ptr_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 425 425 426 rc_pwatt = rc_pwatt * r au0_rcp ! conversion from K.s-1 to PetaWatt426 rc_pwatt = rc_pwatt * rho0_rcp ! conversion from K.s-1 to PetaWatt 427 427 428 428 IF( lk_mpp ) CALL mpp_ini_znl( numout ) ! Define MPI communicator for zonal sum … … 448 448 ! Initialise arrays to zero because diatpr is called before they are first calculated 449 449 ! Note that this means diagnostics will not be exactly correct when model run is restarted. 450 htr_adv(:,:) = 0._wp ; str_adv(:,:) = 0._wp451 htr_ldf(:,:) = 0._wp ; str_ldf(:,:) = 0._wp452 htr_eiv(:,:) = 0._wp ; str_eiv(:,:) = 0._wp450 htr_adv(:,:) = 0._wp ; str_adv(:,:) = 0._wp 451 htr_ldf(:,:) = 0._wp ; str_ldf(:,:) = 0._wp 452 htr_eiv(:,:) = 0._wp ; str_eiv(:,:) = 0._wp 453 453 htr_ove(:,:) = 0._wp ; str_ove(:,:) = 0._wp 454 454 htr_btr(:,:) = 0._wp ; str_btr(:,:) = 0._wp -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diawri.F90
r9652 r9923 169 169 170 170 IF ( iom_use("taubot") ) THEN ! bottom stress 171 zztmp = r au0 * 0.25171 zztmp = rho0 * 0.25 172 172 z2d(:,:) = 0._wp 173 173 DO jj = 2, jpjm1 … … 212 212 IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value 213 213 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 214 z2d(:,:) = r au0 * e1e2t(:,:)214 z2d(:,:) = rho0 * e1e2t(:,:) 215 215 DO jk = 1, jpk 216 216 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) … … 253 253 END DO 254 254 END DO 255 CALL iom_put( "heatc", r au0_rcp * z2d ) ! vertically integrated heat content (J/m2)255 CALL iom_put( "heatc", rho0_rcp * z2d ) ! vertically integrated heat content (J/m2) 256 256 ENDIF 257 257 … … 265 265 END DO 266 266 END DO 267 CALL iom_put( "saltc", r au0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)267 CALL iom_put( "saltc", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2) 268 268 ENDIF 269 269 ! … … 291 291 z2d(:,:) = 0.e0 292 292 DO jk = 1, jpkm1 293 z3d(:,:,jk) = r au0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)293 z3d(:,:,jk) = rho0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk) 294 294 z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 295 295 END DO … … 328 328 z3d(:,:,jpk) = 0.e0 329 329 DO jk = 1, jpkm1 330 z3d(:,:,jk) = r au0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)330 z3d(:,:,jk) = rho0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 331 331 END DO 332 332 CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction … … 369 369 END DO 370 370 CALL lbc_lnk( z2d, 'T', -1. ) 371 CALL iom_put( "tosmint", r au0 * z2d ) ! Vertical integral of temperature371 CALL iom_put( "tosmint", rho0 * z2d ) ! Vertical integral of temperature 372 372 ENDIF 373 373 IF( iom_use("somint") ) THEN … … 381 381 END DO 382 382 CALL lbc_lnk( z2d, 'T', -1. ) 383 CALL iom_put( "somint", r au0 * z2d ) ! Vertical integral of salinity383 CALL iom_put( "somint", rho0 * z2d ) ! Vertical integral of salinity 384 384 ENDIF 385 385 … … 458 458 clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes) 459 459 #if defined key_diainstant 460 zsto = nwrite * r dt460 zsto = nwrite * rn_Dt 461 461 clop = "inst("//TRIM(clop)//")" 462 462 #else 463 zsto =rdt463 zsto = rn_Dt 464 464 clop = "ave("//TRIM(clop)//")" 465 465 #endif 466 zout = nwrite * r dt467 zmax = ( nitend - nit000 + 1 ) * r dt466 zout = nwrite * rn_Dt 467 zmax = ( nitend - nit000 + 1 ) * rn_Dt 468 468 469 469 ! Define indices of the horizontal output zoom and vertical limit storage … … 485 485 486 486 ! Compute julian date from starting date of the run 487 CALL ymds2ju( nyear, nmonth, nday, r dt, zjulian )487 CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian ) 488 488 zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment 489 489 IF(lwp)WRITE(numout,*) … … 507 507 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit 508 508 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 509 & nit000-1, zjulian, r dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )509 & nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 510 510 CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept 511 511 & "m", ipk, gdept_1d, nz_T, "down" ) … … 543 543 CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu 544 544 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 545 & nit000-1, zjulian, r dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )545 & nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 546 546 CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept 547 547 & "m", ipk, gdept_1d, nz_U, "down" ) … … 556 556 CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv 557 557 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 558 & nit000-1, zjulian, r dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )558 & nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 559 559 CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept 560 560 & "m", ipk, gdept_1d, nz_V, "down" ) … … 569 569 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit 570 570 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 571 & nit000-1, zjulian, r dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )571 & nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set ) 572 572 CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw 573 573 & "m", ipk, gdepw_1d, nz_W, "down" ) … … 897 897 clname = cdfile_name 898 898 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 899 zsto = r dt899 zsto = rn_Dt 900 900 clop = "inst(x)" ! no use of the mask value (require less cpu time) 901 zout = r dt902 zmax = ( nitend - nit000 + 1 ) * r dt901 zout = rn_Dt 902 zmax = ( nitend - nit000 + 1 ) * rn_Dt 903 903 904 904 IF(lwp) WRITE(numout,*) … … 912 912 913 913 ! Compute julian date from starting date of the run 914 CALL ymds2ju( nyear, nmonth, nday, r dt, zjulian ) ! time axis914 CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian ) ! time axis 915 915 zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment 916 916 CALL histbeg( clname, jpi, glamt, jpj, gphit, & 917 1, jpi, 1, jpj, nit000-1, zjulian, r dt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit917 1, jpi, 1, jpj, nit000-1, zjulian, rn_Dt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit 918 918 CALL histvert( id_i, "deptht", "Vertical T levels", & ! Vertical grid : gdept 919 919 "m", jpk, gdept_1d, nz_i, "down") -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIU/cool_skin.F90
r9598 r9923 68 68 69 69 70 SUBROUTINE diurnal_sst_coolskin_step(psqflux, pstauflux, psrho, rdt)70 SUBROUTINE diurnal_sst_coolskin_step(psqflux, pstauflux, psrho, pdt) 71 71 !!---------------------------------------------------------------------- 72 72 !! *** ROUTINE diurnal_sst_takaya_step *** … … 82 82 REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux ! Wind stress (kg/ m s^2) 83 83 REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho ! Water density (kg/m^3) 84 REAL(wp), INTENT(IN) :: rdt ! Time-step84 REAL(wp), INTENT(IN) :: pdt ! Time-step (s) 85 85 86 86 ! Local variables -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIU/diurnal_bulk.F90
r9168 r9923 78 78 79 79 80 SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, p _rdt, &80 SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, pdt, & 81 81 & pla, pthick, pcoolthick, pmu, & 82 82 & p_fvel_bkginc, p_hflux_bkginc) … … 98 98 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2) 99 99 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3) 100 REAL(wp) , INTENT(in) :: p _rdt ! time-step100 REAL(wp) , INTENT(in) :: pdt ! time-step (s) 101 101 REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pLa ! Langmuir number 102 102 REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m) … … 167 167 168 168 ! Increment the temperature using the implicit solution 169 x_dsst(:,:) = t_imp( x_dsst(:,:), p _rdt, z_abflux(:,:), z_fvel(:,:), &169 x_dsst(:,:) = t_imp( x_dsst(:,:), pdt, z_abflux(:,:), z_fvel(:,:), & 170 170 & z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:) ) 171 171 ! … … 173 173 174 174 175 FUNCTION t_imp(p_dsst, p _rdt, p_abflux, p_fvel, &175 FUNCTION t_imp(p_dsst, pdt, p_abflux, p_fvel, & 176 176 p_fla, pmu, pthick, prho ) 177 177 … … 182 182 ! Dummy variables 183 183 REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst ! Delta SST 184 REAL(wp), INTENT(IN) :: p _rdt ! Time-step184 REAL(wp), INTENT(IN) :: pdt ! Time-step (s) 185 185 REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux ! Heat forcing 186 186 REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel ! Friction velocity … … 257 257 & ( pthick(ji,jj) * z_stabfunc ) ) 258 258 259 t_imp(ji,jj) = ( p_dsst(ji,jj) + p _rdt * z_term1 ) / &260 ( 1._wp - p _rdt * z_term2 )259 t_imp(ji,jj) = ( p_dsst(ji,jj) + pdt * z_term1 ) / & 260 ( 1._wp - pdt * z_term2 ) 261 261 262 262 END DO 263 263 END DO 264 264 265 END FUNCTION t_imp 266 265 END FUNCTION t_imp 266 267 !!====================================================================== 267 268 END MODULE diurnal_bulk -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIU/step_diu.F90
r9598 r9923 5 5 !!====================================================================== 6 6 !! History : 3.7 ! 2015-11 (J. While) Original code 7 !!---------------------------------------------------------------------- 7 8 8 9 USE diurnal_bulk ! diurnal SST bulk routines (diurnal_sst_takaya routine) … … 27 28 !! Software governed by the CeCILL licence (./LICENSE) 28 29 !!---------------------------------------------------------------------- 29 30 30 CONTAINS 31 31 32 32 SUBROUTINE stp_diurnal( kstp ) 33 INTEGER, INTENT(in) :: kstp ! ocean time-step index34 33 !!---------------------------------------------------------------------- 35 34 !! *** ROUTINE stp_diurnal *** … … 46 45 !! -8- Outputs and diagnostics 47 46 !!---------------------------------------------------------------------- 47 INTEGER, INTENT(in) :: kstp ! ocean time-step index 48 ! 48 49 INTEGER :: jk ! dummy loop indices 49 50 INTEGER :: indic ! error indicator if < 0 … … 51 52 !! --------------------------------------------------------------------- 52 53 53 IF( ln_diurnal_only) THEN54 IF( ln_diurnal_only ) THEN 54 55 indic = 0 ! reset to no error condition 55 56 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) … … 60 61 ENDIF 61 62 62 CALL sbc ( kstp ) ! SeaBoundary Conditions63 CALL sbc( kstp ) ! Sea Surface Boundary Conditions 63 64 ENDIF 64 65 65 ! Cool skin66 66 IF( .NOT.ln_diurnal ) CALL ctl_stop( "stp_diurnal: ln_diurnal not set" ) 67 67 68 68 IF( .NOT. ln_blk ) CALL ctl_stop( "stp_diurnal: diurnal flux processing only implemented for bulk forcing" ) 69 69 70 CALL diurnal_sst_coolskin_step( qns, taum, rhop(:,:,1), rdt) 70 ! ! Cool skin 71 CALL diurnal_sst_coolskin_step( qns, taum, rhop(:,:,1), rn_Dt ) 71 72 72 CALL iom_put( "sst_wl" , x_dsst )! warm layer (write out before update below).73 CALL iom_put( "sst_cs" , x_csdsst )! cool skin73 CALL iom_put( "sst_wl", x_dsst ) ! warm layer (write out before update below). 74 CALL iom_put( "sst_cs", x_csdsst ) ! cool skin 74 75 75 ! Diurnal warm layer model 76 CALL diurnal_sst_takaya_step( kstp, & 77 & qsr, qns, taum, rhop(:,:,1), rdt) 76 ! ! Diurnal warm layer model 77 CALL diurnal_sst_takaya_step( kstp, qsr, qns, taum, rhop(:,:,1), rn_Dt ) 78 78 79 79 IF( ln_diurnal_only ) THEN 80 IF( ln_diaobs ) CALL dia_obs( kstp )! obs-minus-model (assimilation) diagnostics (call after dynamics update)80 IF( ln_diaobs ) CALL dia_obs( kstp ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 81 81 82 82 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 84 84 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 85 85 IF( kstp == nit000 ) CALL iom_close( numror ) ! close input ocean restart file 86 IF( lrst_oce ) CALL rst_write ( kstp )! write output ocean restart file86 IF( lrst_oce ) CALL rst_write( kstp ) ! write output ocean restart file 87 87 88 88 IF( ln_timing .AND. kstp == nit000 ) CALL timing_reset … … 91 91 END SUBROUTINE stp_diurnal 92 92 93 !!====================================================================== 93 94 END MODULE step_diu -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/daymod.F90
r9598 r9923 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_MLF/src/OCE/DOM/dom_oce.F90
r9863 r9923 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 37 LOGICAL , PUBLIC :: ln_1st_euler !: =0 start with forward time step or not (=1) … … 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 tracer (=rn_rdt)60 61 56 ! !!! associated variables 62 57 LOGICAL , PUBLIC :: l_1st_euler !: Euler 1st time-step flag (=T if ln_restart=F or ln_1st_euler=T) 63 REAL(wp), PUBLIC :: r2dt, r1_2dt !: = 2*rdt and 1/(2*rdt) except if l_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_MLF/src/OCE/DOM/domain.F90
r9863 r9923 293 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 … … 414 414 WRITE(numout,*) ' create mesh/mask file ln_meshmask = ', ln_meshmask 415 415 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' [m]' 416 WRITE(numout,*) ' ocean time step rn_ rdt = ', rn_rdt416 WRITE(numout,*) ' ocean time step rn_dt = ', rn_dt 417 417 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 418 418 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs 419 419 ENDIF 420 420 ! 421 ! ! conversion DOCTOR names into model names (this should disappear soon)422 atfp = rn_atfp423 rdt = rn_rdt424 425 421 IF( TRIM(Agrif_CFixed()) == '0' ) THEN 426 422 lrxios = ln_xios_read.AND.ln_rstart -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/domvvl.F90
r9863 r9923 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 ! … … 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 … … 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 … … 607 605 DO ji = 1, jpi 608 606 ze3f = te3t_n(ji,jj,jk) & 609 & + atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) )607 & + rn_atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) ) 610 608 te3t_b(ji,jj,jk) = ze3f 611 609 te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk) … … 1018 1016 WRITE(numout,*) ' rn_rst_e3t = 0.e0' 1019 1017 WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 1020 WRITE(numout,*) ' rn_lf_cutoff = 1 .0/rdt'1018 WRITE(numout,*) ' rn_lf_cutoff = 1/rn_Dt' 1021 1019 ELSE 1022 1020 WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/iscplini.F90
r9598 r9923 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_MLF/src/OCE/DOM/phycst.F90
r9656 r9923 46 46 REAL(wp), PUBLIC :: rt0_ice = 273.05_wp !: melting point of ice [Kelvin] 47 47 #endif 48 REAL(wp), PUBLIC :: r au0 !: volumic mass of reference [kg/m3]49 REAL(wp), PUBLIC :: r1_r au0 !: = 1. / rau0 [m3/kg]48 REAL(wp), PUBLIC :: rho0 !: volumic mass of reference [kg/m3] 49 REAL(wp), PUBLIC :: r1_rho0 !: = 1. / rho0 [m3/kg] 50 50 REAL(wp), PUBLIC :: rcp !: ocean specific heat [J/Kelvin] 51 51 REAL(wp), PUBLIC :: r1_rcp !: = 1. / rcp [Kelvin/J] 52 REAL(wp), PUBLIC :: r au0_rcp !: = rau0 * rcp53 REAL(wp), PUBLIC :: r1_r au0_rcp !: = 1. / ( rau0 * rcp )52 REAL(wp), PUBLIC :: rho0_rcp !: = rho0 * rcp 53 REAL(wp), PUBLIC :: r1_rho0_rcp !: = 1. / ( rho0 * rcp ) 54 54 55 55 REAL(wp), PUBLIC :: rhosn = 330._wp !: volumic mass of snow [kg/m3] -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/restart.F90
r9863 r9923 148 148 IF( lwxios ) CALL iom_swap( cwxios_context ) 149 149 150 CALL iom_rstput( kt, nitrst, numrow, 'rdt' , rn_ rdt, ldxios = lwxios ) ! dynamics time step150 CALL iom_rstput( kt, nitrst, numrow, 'rdt' , rn_Dt , ldxios = lwxios ) ! dynamics time step 151 151 ! 152 152 IF( .NOT. ln_diurnal_only ) THEN … … 261 261 IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 ) THEN 262 262 CALL iom_get( numror, 'rdt', zrdt, ldxios = lrxios ) 263 IF( zrdt /= rn_ rdt ) THEN263 IF( zrdt /= rn_Dt ) THEN 264 264 IF(lwp) WRITE( numout,*) 265 265 IF(lwp) WRITE( numout,*) 'rst_read: rdt not equal to the read one' … … 273 273 IF( ln_diurnal ) CALL iom_get( numror, jpdom_autoglo, 'Dsst' , x_dsst, ldxios = lrxios ) 274 274 IF ( ln_diurnal_only ) THEN 275 IF(lwp) WRITE( numout,*) 'rst_read: ln_diurnal_only set, setting rhop=r au0'276 rhop = r au0275 IF(lwp) WRITE( numout,*) 'rst_read: ln_diurnal_only set, setting rhop=rho0' 276 rhop = rho0 277 277 CALL iom_get( numror, jpdom_autoglo, 'tn' , w3d, ldxios = lrxios ) 278 278 tsn(:,:,1,jp_tem) = w3d(:,:,1) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/divhor.F90
r9598 r9923 63 63 ! 64 64 INTEGER :: ji, jj, jk ! dummy loop indices 65 REAL(wp) :: z raur, zdep! local scalars65 REAL(wp) :: zdep ! local scalars 66 66 !!---------------------------------------------------------------------- 67 67 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynnxt.F90
r9863 r9923 83 83 !! * Apply the time filter applied and swap of the dynamics 84 84 !! arrays to start the next time step: 85 !! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]85 !! (ub,vb) = (un,vn) + rn_atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 86 86 !! (un,vn) = (ua,va). 87 87 !! Note that with flux form advection and non linear free surface, … … 156 156 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends 157 157 IF( ln_dynadv_vec ) THEN 158 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_ 2dt159 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_ 2dt158 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt 159 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt 160 160 ELSE 161 zua(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_ 2dt162 zva(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_ 2dt161 zua(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt 162 zva(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt 163 163 ENDIF 164 164 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin filter … … 192 192 DO jj = 1, jpj 193 193 DO ji = 1, jpi 194 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )195 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )194 zuf = un(ji,jj,jk) + rn_atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 195 zvf = vn(ji,jj,jk) + rn_atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 196 196 ! 197 197 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 207 207 ! Before scale factor at t-points (used as a now filtered scale factor until the swap) 208 208 DO jk = 1, jpkm1 209 e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )209 e3t_b(:,:,jk) = e3t_n(:,:,jk) + rn_atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 210 210 END DO 211 211 ! Add volume filter correction: compatibility with tracer advection scheme 212 212 ! => time filter + conservation correction (only at the first level) 213 zcoef = atfp * rdt * r1_rau0213 zcoef = rn_atfp * rn_Dt * r1_rho0 214 214 215 215 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) … … 252 252 DO jj = 1, jpj 253 253 DO ji = 1, jpi 254 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )255 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )254 zuf = un(ji,jj,jk) + rn_atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 255 zvf = vn(ji,jj,jk) + rn_atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 256 256 ! 257 257 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 279 279 zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) 280 280 ! 281 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)282 zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)281 zuf = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 282 zvf = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 283 283 ! 284 284 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 349 349 ENDIF 350 350 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 351 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * r1_ 2dt352 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * r1_ 2dt351 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * r1_Dt 352 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * r1_Dt 353 353 CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 354 354 ENDIF -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg.F90
r9863 r9923 66 66 !! ln_apr_dyn=T : the atmospheric pressure forcing is applied 67 67 !! as the gradient of the inverse barometer ssh: 68 !! apgu = - 1/r au0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]69 !! apgv = - 1/r au0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]68 !! apgu = - 1/rho0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb] 69 !! apgv = - 1/rho0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb] 70 70 !! Note that as all external forcing a time averaging over a two rdt 71 71 !! period is used to prevent the divergence of odd and even time step. … … 74 74 ! 75 75 INTEGER :: ji, jj, jk ! dummy loop indices 76 REAL(wp) :: zg_2, zintp, zg rau0r, zld ! local scalars76 REAL(wp) :: zg_2, zintp, zg_rho0, zld ! local scalars 77 77 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice 78 78 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv … … 110 110 ENDIF 111 111 ! 112 ! !== tide potential forcing term ==! 113 IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. ln_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case 114 ! 115 CALL upd_tide( kt ) ! update tide potential 116 ! 117 DO jj = 2, jpjm1 ! add tide potential forcing 118 DO ji = fs_2, fs_jpim1 ! vector opt. 119 spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 120 spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 121 END DO 122 END DO 123 ! 124 IF (ln_scal_load) THEN 125 zld = rn_scal_load * grav 126 DO jj = 2, jpjm1 ! add scalar approximation for load potential 127 DO ji = fs_2, fs_jpim1 ! vector opt. 128 spgu(ji,jj) = spgu(ji,jj) + zld * ( sshn(ji+1,jj) - sshn(ji,jj) ) * r1_e1u(ji,jj) 129 spgv(ji,jj) = spgv(ji,jj) + zld * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_e2v(ji,jj) 130 END DO 131 END DO 112 IF( .NOT.ln_dynspg_ts ) THEN 113 ! !== tide potential forcing term ==! 114 IF( ln_tide_pot .AND. ln_tide ) THEN ! N.B. added directly at sub-time-step in ts-case 115 ! 116 CALL upd_tide( kt ) ! update tide potential 117 ! 118 IF ( ln_scal_load ) THEN 119 zld = rn_load * grav 120 DO jj = 2, jpjm1 ! add tide potential + scalar approximation of load potential 121 DO ji = fs_2, fs_jpim1 ! vector opt. 122 spgu(ji,jj) = spgu(ji,jj) + ( grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) & 123 & + zld * ( sshn (ji+1,jj) - sshn (ji,jj) ) ) * r1_e1u(ji,jj) 124 spgv(ji,jj) = spgv(ji,jj) + ( grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) & 125 & + zld * ( sshn (ji,jj+1) - sshn (ji,jj) ) ) * r1_e2v(ji,jj) 126 END DO 127 END DO 128 ELSE 129 DO jj = 2, jpjm1 ! add tide potential 130 DO ji = fs_2, fs_jpim1 ! vector opt. 131 spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 132 spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 133 END DO 134 END DO 135 ENDIF 132 136 ENDIF 133 137 ENDIF … … 136 140 ALLOCATE( zpice(jpi,jpj) ) 137 141 zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc ) 138 zg rau0r = - grav * r1_rau0139 zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zg rau0r142 zg_rho0 = - grav * r1_rho0 143 zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zg_rho0 140 144 DO jj = 2, jpjm1 141 145 DO ji = fs_2, fs_jpim1 ! vector opt. … … 191 195 NAMELIST/namdyn_spg/ ln_dynspg_exp , ln_dynspg_ts, & 192 196 & ln_bt_fw, ln_bt_av , ln_bt_auto , & 193 & nn_ baro, rn_bt_cmax, nn_bt_flt, rn_bt_alpha197 & nn_e , rn_bt_cmax, nn_bt_flt, rn_bt_alpha 194 198 !!---------------------------------------------------------------------- 195 199 ! … … 227 231 WRITE(numout,*) 228 232 IF( nspg == np_EXP ) WRITE(numout,*) ' ==>>> explicit free surface' 229 IF( nspg == np_TS ) WRITE(numout,*) ' ==>>> free surface with time splitting scheme'233 IF( nspg == np_TS ) WRITE(numout,*) ' ==>>> split-explicit free surface' 230 234 IF( nspg == np_NO ) WRITE(numout,*) ' ==>>> No surface surface pressure gradient trend in momentum Eqs.' 231 235 ENDIF 232 236 ! 233 237 IF( nspg == np_TS ) THEN ! split-explicit scheme initialisation 234 CALL dyn_spg_ts_init ! do it first: set nn_ baroused to allocate some arrays later on238 CALL dyn_spg_ts_init ! do it first: set nn_e used to allocate some arrays later on 235 239 ENDIF 236 240 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg_exp.F90
r9598 r9923 49 49 !! momentum trend the surface pressure gradient : 50 50 !! (ua,va) = (ua,va) + (spgu,spgv) 51 !! where spgu = -1/r au0 d/dx(ps) = -g/e1u di( sshn )52 !! spgv = -1/r au0 d/dy(ps) = -g/e2v dj( sshn )51 !! where spgu = -1/rho0 d/dx(ps) = -g/e1u di( sshn ) 52 !! spgv = -1/rho0 d/dy(ps) = -g/e2v dj( sshn ) 53 53 !! 54 54 !! ** Action : (ua,va) trend of horizontal velocity increased by -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg_ts.F90
r9863 r9923 69 69 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 70 70 ! 71 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_ baro <= 2.5 nn_baro72 REAL(wp),SAVE :: r dtbt ! Barotropic timestep71 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5*nn_e 72 REAL(wp),SAVE :: rDt_e ! external mode time-step 73 73 ! 74 74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields … … 81 81 REAL(wp) :: r1_4 = 0.25_wp ! 82 82 REAL(wp) :: r1_2 = 0.5_wp ! 83 84 REAL(wp) :: r1_2dt_b, r2dt_bf ! local scalars85 83 86 84 !! * Substitutions … … 101 99 ierr(:) = 0 102 100 ! 103 ALLOCATE( wgtbtp1(3*nn_ baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )101 ALLOCATE( wgtbtp1(3*nn_e), wgtbtp2(3*nn_e), zwz(jpi,jpj), STAT=ierr(1) ) 104 102 ! 105 103 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & … … 150 148 INTEGER :: ikbu, iktu, noffset ! local integers 151 149 INTEGER :: ikbv, iktv ! - - 150 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point 152 151 REAL(wp) :: zx1, zx2, zu_spg, zhura, z1_hu ! - - 153 152 REAL(wp) :: zy1, zy2, zv_spg, zhvra, z1_hv ! - - 154 153 REAL(wp) :: za0, za1, za2, za3 ! - - 155 154 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 155 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 156 REAL(wp) :: zload 156 157 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 157 158 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc … … 161 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 162 163 ! 163 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 164 165 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point 164 166 165 167 166 REAL(wp) :: zepsilon, zgamma ! - - … … 179 178 zwdramp = r_rn_wdmin1 ! simplest ramp 180 179 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 181 ! ! reciprocal of baroclinic time step182 IF( l_1st_euler ) THEN ; r2dt_bf = rdt183 ELSE ; r2dt_bf = 2.0_wp * rdt184 ENDIF185 r1_2dt_b = 1.0_wp / r2dt_bf186 180 ! 187 181 ll_init = ln_bt_av ! if no time averaging, then no specific restart 188 182 ll_fw_start = .FALSE. 189 183 ! ! time offset in steps for bdy data update 190 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_ baro184 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e 191 185 ELSE ; noffset = 0 192 186 ENDIF … … 459 453 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 460 454 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 461 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)455 zcpx(ji,jj) = MAX( 0._wp , MIN( zcpx(ji,jj) , 1._wp ) ) 462 456 ELSE 463 457 zcpx(ji,jj) = 0._wp … … 536 530 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 537 531 IF( ln_wd_il ) THEN 538 zztmp = -1._wp / r dtbt532 zztmp = -1._wp / rDt_e 539 533 DO jj = 2, jpjm1 540 534 DO ji = fs_2, fs_jpim1 ! vector opt. … … 587 581 DO jj = 2, jpjm1 588 582 DO ji = fs_2, fs_jpim1 ! vector opt. 589 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_r au0 * utau(ji,jj) * r1_hu_n(ji,jj)590 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_r au0 * vtau(ji,jj) * r1_hv_n(ji,jj)583 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu_n(ji,jj) 584 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv_n(ji,jj) 591 585 END DO 592 586 END DO 593 587 ELSE 594 zztmp = r1_r au0 * r1_2588 zztmp = r1_rho0 * r1_2 595 589 DO jj = 2, jpjm1 596 590 DO ji = fs_2, fs_jpim1 ! vector opt. … … 629 623 ! ! Surface net water flux and rivers 630 624 IF (ln_bt_fw) THEN 631 zssh_frc(:,:) = r1_r au0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )625 zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 632 626 ELSE 633 zztmp = r1_r au0 * r1_2627 zztmp = r1_rho0 * r1_2 634 628 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 635 629 & + fwfisf(:,:) + fwfisf_b(:,:) ) … … 818 812 ENDIF 819 813 #endif 820 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, r dtbt)814 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rDt_e) 821 815 822 816 IF ( ln_wd_dl ) THEN … … 864 858 END DO 865 859 END DO 866 ssha_e(:,:) = ( sshn_e(:,:) - r dtbt* ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)860 ssha_e(:,:) = ( sshn_e(:,:) - rDt_e * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 867 861 868 862 CALL lbc_lnk( ssha_e, 'T', 1._wp ) … … 1068 1062 ENDIF 1069 1063 ! 1070 ! Surface pressure trend: 1064 ! Surface pressure trend 1065 IF( ln_scal_load ) THEN ; zload = 1._wp 1066 ELSE ; zload = 1._wp - rn_load 1067 ENDIF 1071 1068 IF( ln_wd_il ) THEN 1072 1069 DO jj = 2, jpjm1 … … 1075 1072 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1076 1073 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1077 zwx(ji,jj) = (1._wp - rn_scal_load)* zu_spg * zcpx(ji,jj)1078 zwy(ji,jj) = (1._wp - rn_scal_load)* zv_spg * zcpy(ji,jj)1074 zwx(ji,jj) = zload * zu_spg * zcpx(ji,jj) 1075 zwy(ji,jj) = zload * zv_spg * zcpy(ji,jj) 1079 1076 END DO 1080 1077 END DO … … 1085 1082 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1086 1083 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1087 zwx(ji,jj) = (1._wp - rn_scal_load)* zu_spg1088 zwy(ji,jj) = (1._wp - rn_scal_load)* zv_spg1084 zwx(ji,jj) = zload * zu_spg 1085 zwy(ji,jj) = zload * zv_spg 1089 1086 END DO 1090 1087 END DO … … 1097 1094 DO ji = fs_2, fs_jpim1 ! vector opt. 1098 1095 ua_e(ji,jj) = ( un_e(ji,jj) & 1099 & + r dtbt* ( zwx(ji,jj) &1096 & + rDt_e * ( zwx(ji,jj) & 1100 1097 & + zu_trd(ji,jj) & 1101 1098 & + zu_frc(ji,jj) ) & … … 1103 1100 1104 1101 va_e(ji,jj) = ( vn_e(ji,jj) & 1105 & + r dtbt* ( zwy(ji,jj) &1102 & + rDt_e * ( zwy(ji,jj) & 1106 1103 & + zv_trd(ji,jj) & 1107 1104 & + zv_frc(ji,jj) ) & … … 1110 1107 !jth implicit bottom friction: 1111 1108 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 1112 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - r dtbt* zCdU_u(ji,jj) * hur_e(ji,jj))1113 va_e(ji,jj) = va_e(ji,jj) /(1.0 - r dtbt* zCdU_v(ji,jj) * hvr_e(ji,jj))1109 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 1110 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) 1114 1111 ENDIF 1115 1112 … … 1128 1125 1129 1126 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 1130 & + r dtbt* ( zhust_e(ji,jj) * zwx(ji,jj) &1127 & + rDt_e * ( zhust_e(ji,jj) * zwx(ji,jj) & 1131 1128 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 1132 1129 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & … … 1134 1131 1135 1132 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 1136 & + r dtbt* ( zhvst_e(ji,jj) * zwy(ji,jj) &1133 & + rDt_e * ( zhvst_e(ji,jj) * zwy(ji,jj) & 1137 1134 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 1138 1135 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & … … 1202 1199 zwy(:,:) = vn_adv(:,:) 1203 1200 IF( .NOT.l_1st_euler ) THEN 1204 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )1205 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )1201 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - rn_atfp * un_bf(:,:) ) 1202 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - rn_atfp * vn_bf(:,:) ) 1206 1203 ! 1207 1204 ! Update corrective fluxes for next time step: 1208 un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))1209 vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))1205 un_bf(:,:) = rn_atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 1206 vn_bf(:,:) = rn_atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 1210 1207 ELSE 1211 1208 un_bf(:,:) = 0._wp … … 1222 1219 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1223 1220 DO jk=1,jpkm1 1224 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_ 2dt_b1225 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_ 2dt_b1221 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_Dt 1222 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_Dt 1226 1223 END DO 1227 1224 ELSE … … 1229 1226 DO jj = 1, jpjm1 1230 1227 DO ji = 1, jpim1 ! NO Vector Opt. 1231 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 1232 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1233 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1234 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 1235 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1236 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1228 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1229 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1230 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1231 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1237 1232 END DO 1238 1233 END DO … … 1240 1235 ! 1241 1236 DO jk=1,jpkm1 1242 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_ 2dt_b1243 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_ 2dt_b1237 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_Dt 1238 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_Dt 1244 1239 END DO 1245 1240 ! Save barotropic velocities not transport: … … 1306 1301 LOGICAL , INTENT(in ) :: ll_fw ! forward time splitting =.true. 1307 1302 INTEGER , INTENT(inout) :: jpit ! cycle length 1308 REAL(wp), DIMENSION(3*nn_ baro), INTENT(inout) :: zwgt1, zwgt2 ! Primary & Secondary weights1303 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, zwgt2 ! Primary & Secondary weights 1309 1304 ! 1310 1305 INTEGER :: jic, jn, ji ! temporary integers … … 1316 1311 ! 1317 1312 ! Set time index when averaged value is requested 1318 IF ( ll_fw ) THEN ; jic = nn_ baro1319 ELSE ; jic = 2 * nn_ baro1313 IF ( ll_fw ) THEN ; jic = nn_e 1314 ELSE ; jic = 2 * nn_e 1320 1315 ENDIF 1321 1316 … … 1323 1318 ! 1324 1319 IF (ll_av) THEN !* Define simple boxcar window for primary weights 1325 ! ! (width = nn_ baro, centered around jic)1320 ! ! (width = nn_e, centered around jic) 1326 1321 SELECT CASE ( nn_bt_flt ) 1327 1322 CASE( 0 ) ! No averaging … … 1329 1324 jpit = jic 1330 1325 ! 1331 CASE( 1 ) ! Boxcar, width = nn_ baro1332 DO jn = 1, 3*nn_ baro1333 za1 = ABS(float(jn-jic))/float(nn_ baro)1326 CASE( 1 ) ! Boxcar, width = nn_e 1327 DO jn = 1, 3*nn_e 1328 za1 = ABS(float(jn-jic))/float(nn_e) 1334 1329 IF ( za1 < 0.5_wp ) THEN 1335 1330 zwgt1(jn) = 1._wp … … 1338 1333 END DO 1339 1334 ! 1340 CASE( 2 ) ! Boxcar, width = 2 * nn_ baro1341 DO jn = 1, 3*nn_ baro1342 za1 = ABS(float(jn-jic))/float(nn_ baro)1335 CASE( 2 ) ! Boxcar, width = 2 * nn_e 1336 DO jn = 1, 3*nn_e 1337 za1 = ABS(float(jn-jic))/float(nn_e) 1343 1338 IF ( za1 < 1._wp ) THEN 1344 1339 zwgt1(jn) = 1._wp … … 1474 1469 1475 1470 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1476 IF( ln_bt_auto ) nn_ baro = CEILING( rdt / rn_bt_cmax * zcmax)1471 IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax) 1477 1472 1478 r dtbt = rdt / REAL( nn_baro, wp )1479 zcmax = zcmax * r dtbt1473 rDt_e = rn_Dt / REAL( nn_e , wp ) 1474 zcmax = zcmax * rDt_e 1480 1475 ! Print results 1481 1476 IF(lwp) WRITE(numout,*) … … 1483 1478 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 1484 1479 IF( ln_bt_auto ) THEN 1485 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_ baro'1480 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e ' 1486 1481 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1487 1482 ELSE 1488 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_ baro in namelist nn_baro = ', nn_baro1483 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e 1489 1484 ENDIF 1490 1485 1491 1486 IF(ln_bt_av) THEN 1492 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_ barotime steps is on '1487 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on ' 1493 1488 ELSE 1494 1489 IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables ' … … 1510 1505 SELECT CASE ( nn_bt_flt ) 1511 1506 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1512 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_ baro'1513 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_ baro'1507 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e' 1508 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1514 1509 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' ) 1515 1510 END SELECT 1516 1511 ! 1517 1512 IF(lwp) WRITE(numout,*) ' ' 1518 IF(lwp) WRITE(numout,*) ' nn_ baro = ', nn_baro1519 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt1520 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax1513 IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e 1514 IF(lwp) WRITE(numout,*) ' external mode time step is : rDt_e', rDt_e, ' [s]' 1515 IF(lwp) WRITE(numout,*) ' Maximum Courant number is : ', zcmax 1521 1516 ! 1522 1517 IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha … … 1529 1524 ENDIF 1530 1525 IF( zcmax>0.9_wp ) THEN 1531 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_ baro!' )1526 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1532 1527 ENDIF 1533 1528 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynzdf.F90
r9863 r9923 90 90 ENDIF 91 91 ! 92 z2dt_2 = r 2dt * 0.5_wp !* =rdt except in 1st Euler time step where it is equal to rdt/292 z2dt_2 = rDt * 0.5_wp !* =rn_Dt except in 1st Euler time step where it is equal to rn_Dt/2 93 93 ! 94 94 ! … … 108 108 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 109 109 DO jk = 1, jpkm1 110 ua(:,:,jk) = ( ub(:,:,jk) + r 2dt * ua(:,:,jk) ) * umask(:,:,jk)111 va(:,:,jk) = ( vb(:,:,jk) + r 2dt * va(:,:,jk) ) * vmask(:,:,jk)110 ua(:,:,jk) = ( ub(:,:,jk) + rDt * ua(:,:,jk) ) * umask(:,:,jk) 111 va(:,:,jk) = ( vb(:,:,jk) + rDt * va(:,:,jk) ) * vmask(:,:,jk) 112 112 END DO 113 113 ELSE ! applied on thickness weighted velocity 114 114 DO jk = 1, jpkm1 115 ua(:,:,jk) = ( e3u_b(:,:,jk)*ub(:,:,jk) + r 2dt * e3u_n(:,:,jk)*ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)116 va(:,:,jk) = ( e3v_b(:,:,jk)*vb(:,:,jk) + r 2dt * e3v_n(:,:,jk)*va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk)115 ua(:,:,jk) = ( e3u_b(:,:,jk)*ub(:,:,jk) + rDt * e3u_n(:,:,jk)*ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 116 va(:,:,jk) = ( e3v_b(:,:,jk)*vb(:,:,jk) + rDt * e3v_n(:,:,jk)*va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 117 117 END DO 118 118 ENDIF … … 160 160 DO ji = fs_2, fs_jpim1 ! vector opt. 161 161 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 162 zzwi = - r 2dt * ( 0.5 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) + akzu(ji,jj,jk ) ) &162 zzwi = - rDt * ( 0.5 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) + akzu(ji,jj,jk ) ) & 163 163 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 164 zzws = - r 2dt * ( 0.5 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) ) &164 zzws = - rDt * ( 0.5 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) ) & 165 165 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 166 166 zwi(ji,jj,jk) = zzwi … … 244 244 DO ji = fs_2, fs_jpim1 ! vector opt. 245 245 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 246 ua(ji,jj,1) = ua(ji,jj,1) + z2dt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( ze3ua * r au0 ) * umask(ji,jj,1)246 ua(ji,jj,1) = ua(ji,jj,1) + z2dt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( ze3ua * rho0 ) * umask(ji,jj,1) 247 247 END DO 248 248 END DO … … 277 277 DO ji = fs_2, fs_jpim1 ! vector opt. 278 278 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 279 zzwi = - r 2dt * ( 0.5 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) + akzv(ji,jj,jk ) ) &279 zzwi = - rDt * ( 0.5 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) + akzv(ji,jj,jk ) ) & 280 280 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 281 zzws = - r 2dt * ( 0.5 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) ) &281 zzws = - rDt * ( 0.5 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) ) & 282 282 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 283 283 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) … … 359 359 DO ji = fs_2, fs_jpim1 ! vector opt. 360 360 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 361 va(ji,jj,1) = va(ji,jj,1) + z2dt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( ze3va * r au0 ) * vmask(ji,jj,1)361 va(ji,jj,1) = va(ji,jj,1) + z2dt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( ze3va * rho0 ) * vmask(ji,jj,1) 362 362 END DO 363 363 END DO … … 385 385 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 386 386 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 387 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_ 2dt - ztrdu(:,:,:)388 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_ 2dt - ztrdv(:,:,:)387 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt - ztrdu(:,:,:) 388 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt - ztrdv(:,:,:) 389 389 ELSE ! applied on thickness weighted velocity 390 ztrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_ 2dt - ztrdu(:,:,:)391 ztrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_ 2dt - ztrdv(:,:,:)390 ztrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt - ztrdu(:,:,:) 391 ztrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt - ztrdv(:,:,:) 392 392 ENDIF 393 393 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) … … 440 440 ! 441 441 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 442 ptrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_ 2dt - ptrdu(:,:,:)443 ptrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_ 2dt - ptrdv(:,:,:)442 ptrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt - ptrdu(:,:,:) 443 ptrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt - ptrdv(:,:,:) 444 444 ELSE ! applied on thickness weighted velocity 445 ptrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_ 2dt - ptrdu(:,:,:)446 ptrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_ 2dt - ptrdv(:,:,:)445 ptrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt - ptrdu(:,:,:) 446 ptrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt - ptrdv(:,:,:) 447 447 ENDIF 448 448 CALL trd_dyn( ptrdu, ptrdv, jpdyn_zdf, kt ) … … 494 494 END DO 495 495 END DO 496 zzz= - 0.5_wp* r au0 ! caution sign minus here496 zzz= - 0.5_wp* rho0 ! caution sign minus here 497 497 z2d(:,:) = zzz * z2d(:,:) 498 498 CALL lbc_lnk( z2d,'T', 1. ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/sshwzv.F90
r9863 r9923 69 69 ! 70 70 INTEGER :: jk ! dummy loop indice 71 REAL(wp) :: z1_2r au0 ! local scalars71 REAL(wp) :: z1_2rho0 ! local scalars 72 72 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 73 73 !!---------------------------------------------------------------------- … … 81 81 ENDIF 82 82 ! 83 z1_2r au0 = 0.5_wp * r1_rau083 z1_2rho0 = 0.5_wp * r1_rho0 84 84 85 85 ! !------------------------------! … … 87 87 ! !------------------------------! 88 88 89 IF(ln_wd_il) CALL wad_lmt( sshb, z1_2r au0 * (emp_b(:,:) + emp(:,:)), r2dt )89 IF(ln_wd_il) CALL wad_lmt( sshb, z1_2rho0 * (emp_b(:,:) + emp(:,:)), rDt ) 90 90 91 91 CALL div_hor( kt ) ! Horizontal divergence … … 99 99 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 100 100 ! 101 ssha(:,:) = ( sshb(:,:) - r 2dt * ( z1_2rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)101 ssha(:,:) = ( sshb(:,:) - rDt * ( z1_2rho0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 102 102 ! 103 103 #if defined key_agrif … … 174 174 ! computation of w 175 175 wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) & 176 & + r1_ 2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )) * tmask(:,:,jk)176 & + r1_Dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk) 177 177 END DO 178 178 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 … … 182 182 ! computation of w 183 183 wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) & 184 & + r1_ 2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)184 & + r1_Dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk) 185 185 END DO 186 186 ENDIF … … 194 194 #if defined key_agrif 195 195 IF( .NOT. AGRIF_Root() ) THEN 196 IF ( (nbondi == 1).OR.(nbondi == 2)) wn(nlci-1 , : ,:) = 0.e0! east197 IF ( (nbondi == -1).OR.(nbondi == 2)) wn(2 , : ,:) = 0.e0! west198 IF ( (nbondj == 1).OR.(nbondj == 2)) wn(: ,nlcj-1 ,:) = 0.e0! north199 IF ( (nbondj == -1).OR.(nbondj == 2)) wn(: ,2 ,:) = 0.e0! south196 IF ( nbondi == 1 .OR. nbondi == 2 ) wn(nlci-1 , : ,:) = 0._wp ! east 197 IF ( nbondi == -1 .OR. nbondi == 2 ) wn( 2 , : ,:) = 0._wp ! west 198 IF ( nbondj == 1 .OR. nbondj == 2 ) wn( : ,nlcj-1 ,:) = 0._wp ! north 199 IF ( nbondj == -1 .OR. nbondj == 2 ) wn( : , 2 ,:) = 0._wp ! south 200 200 ENDIF 201 201 #endif … … 216 216 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 217 217 !! from the filter, see Leclair and Madec 2010) and swap : 218 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha )219 !! - atfp * rdt * ( emp_b - emp ) / rau0218 !! sshn = ssha + rn_atfp * ( sshb -2 sshn + ssha ) 219 !! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0 220 220 !! sshn = ssha 221 221 !! … … 245 245 ! 246 246 ! ! before <-- now filtered 247 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )247 sshb(:,:) = sshn(:,:) + rn_atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 248 248 IF( .NOT.ln_linssh ) THEN ! before <-- with forcing removed 249 zcoef = atfp * rdt * r1_rau0249 zcoef = rn_atfp * rn_Dt * r1_rho0 250 250 sshb(:,:) = sshb(:,:) - zcoef * ( emp_b(:,:) - emp (:,:) & 251 251 & - rnf_b(:,:) + rnf (:,:) & -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/FLO/flo4rk.F90
r9598 r9923 131 131 ! computation of Runge-Kutta factor 132 132 DO jfl = 1, jpnfl 133 zrkxfl(jfl,jind) = r dt*zufl(jfl)134 zrkyfl(jfl,jind) = r dt*zvfl(jfl)135 zrkzfl(jfl,jind) = r dt*zwfl(jfl)133 zrkxfl(jfl,jind) = rn_Dt*zufl(jfl) 134 zrkyfl(jfl,jind) = rn_Dt*zvfl(jfl) 135 zrkzfl(jfl,jind) = rn_Dt*zwfl(jfl) 136 136 END DO 137 137 IF( jind /= 4 ) THEN -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/FLO/floblk.F90
r9598 r9923 234 234 ! test to know if the "age" of the float is not bigger than the 235 235 ! time step 236 IF( zagenewfl(jfl) > r dt ) THEN237 zttfl(jfl) = (r dt-zagefl(jfl)) / zvol238 zagenewfl(jfl) = r dt236 IF( zagenewfl(jfl) > rn_Dt ) THEN 237 zttfl(jfl) = (rn_Dt-zagefl(jfl)) / zvol 238 zagenewfl(jfl) = rn_Dt 239 239 ENDIF 240 240 … … 341 341 ifin = 1 342 342 DO jfl = 1, jpnfl 343 IF( zagefl(jfl) < r dt ) ifin = 0343 IF( zagefl(jfl) < rn_Dt ) ifin = 0 344 344 tpifl(jfl) = zgifl(jfl) + 0.5 345 345 tpjfl(jfl) = zgjfl(jfl) + 0.5 … … 348 348 ifin = 1 349 349 DO jfl = 1, jpnfl 350 IF( zagefl(jfl) < r dt ) ifin = 0350 IF( zagefl(jfl) < rn_Dt ) ifin = 0 351 351 tpifl(jfl) = zgifl(jfl) + 0.5 352 352 tpjfl(jfl) = zgjfl(jfl) + 0.5 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/FLO/flowri.F90
r9598 r9923 125 125 ztem(jfl) = tsn(iafloc,ibfloc,icfl,jp_tem) 126 126 zsal (jfl) = tsn(iafloc,ibfloc,icfl,jp_sal) 127 zrho (jfl) = (rhd(iafloc,ibfloc,icfl)+1)*r au0127 zrho (jfl) = (rhd(iafloc,ibfloc,icfl)+1)*rho0 128 128 129 129 ENDIF … … 145 145 ztem(jfl) = tsn(iafloc,ibfloc,icfl,jp_tem) 146 146 zsal(jfl) = tsn(iafloc,ibfloc,icfl,jp_sal) 147 zrho(jfl) = (rhd(iafloc,ibfloc,icfl)+1)*r au0147 zrho(jfl) = (rhd(iafloc,ibfloc,icfl)+1)*rho0 148 148 149 149 ENDIF … … 248 248 !------------------------------- 249 249 irec = INT( (kt-nn_it000+1)/nn_writefl ) +1 250 ztime = ( kt-nn_it000 + 1 ) * r dt250 ztime = ( kt-nn_it000 + 1 ) * rn_Dt 251 251 252 252 CALL flioputv( numflo , 'time_counter', ztime , start=(/irec/) ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/ICB/icbini.F90
r9598 r9923 58 58 !! - setup either test icebergs or calving file 59 59 !!---------------------------------------------------------------------- 60 REAL(wp), INTENT(in) :: pdt ! iceberg time-step (r dt*nn_fsbc)60 REAL(wp), INTENT(in) :: pdt ! iceberg time-step (rn_Dt*nn_fsbc) 61 61 INTEGER , INTENT(in) :: kt ! time step number 62 62 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/ICB/icbtrj.F90
r9598 r9923 69 69 !!---------------------------------------------------------------------- 70 70 71 !!gm we could probably use the daymod calculation here.... 72 !! ===>>> TO BE checked by someone 73 71 74 ! compute initial time step date 72 75 CALL ju2ymds( fjulday, iyear, imonth, iday, zsec ) … … 74 77 75 78 ! compute end time step date 76 zfjulday = fjulday + r dt / rday * REAL( nitend - nit000 + 1 , wp)79 zfjulday = fjulday + rn_Dt / rday * REAL( nitend - nit000 + 1 , wp) 77 80 IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday ) zfjulday = REAL(NINT(zfjulday),wp) ! avoid truncation error 78 81 CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/IOM/iom.F90
r9802 r9923 239 239 ! 240 240 ! end file definition 241 dtime%second = r dt241 dtime%second = rn_Dt 242 242 CALL xios_set_timestep( dtime ) 243 243 CALL xios_close_context_definition() … … 2358 2358 idx = INDEX(clname,'@startdate@') + INDEX(clname,'@STARTDATE@') 2359 2359 DO WHILE ( idx /= 0 ) 2360 cldate = iom_sdate( fjulday - r dt / rday )2360 cldate = iom_sdate( fjulday - rn_Dt / rday ) 2361 2361 clname = clname(1:idx-1)//TRIM(cldate)//clname(idx+11:LEN_TRIM(clname)) 2362 2362 idx = INDEX(clname,'@startdate@') + INDEX(clname,'@STARTDATE@') … … 2365 2365 idx = INDEX(clname,'@startdatefull@') + INDEX(clname,'@STARTDATEFULL@') 2366 2366 DO WHILE ( idx /= 0 ) 2367 cldate = iom_sdate( fjulday - r dt / rday, ldfull = .TRUE. )2367 cldate = iom_sdate( fjulday - rn_Dt / rday, ldfull = .TRUE. ) 2368 2368 clname = clname(1:idx-1)//TRIM(cldate)//clname(idx+15:LEN_TRIM(clname)) 2369 2369 idx = INDEX(clname,'@startdatefull@') + INDEX(clname,'@STARTDATEFULL@') … … 2372 2372 idx = INDEX(clname,'@enddate@') + INDEX(clname,'@ENDDATE@') 2373 2373 DO WHILE ( idx /= 0 ) 2374 cldate = iom_sdate( fjulday + r dt / rday * REAL( nitend - nit000, wp ), ld24 = .TRUE. )2374 cldate = iom_sdate( fjulday + rn_Dt / rday * REAL( nitend - nit000, wp ), ld24 = .TRUE. ) 2375 2375 clname = clname(1:idx-1)//TRIM(cldate)//clname(idx+9:LEN_TRIM(clname)) 2376 2376 idx = INDEX(clname,'@enddate@') + INDEX(clname,'@ENDDATE@') … … 2379 2379 idx = INDEX(clname,'@enddatefull@') + INDEX(clname,'@ENDDATEFULL@') 2380 2380 DO WHILE ( idx /= 0 ) 2381 cldate = iom_sdate( fjulday + r dt / rday * REAL( nitend - nit000, wp ), ld24 = .TRUE., ldfull = .TRUE. )2381 cldate = iom_sdate( fjulday + rn_Dt / rday * REAL( nitend - nit000, wp ), ld24 = .TRUE., ldfull = .TRUE. ) 2382 2382 clname = clname(1:idx-1)//TRIM(cldate)//clname(idx+13:LEN_TRIM(clname)) 2383 2383 idx = INDEX(clname,'@enddatefull@') + INDEX(clname,'@ENDDATEFULL@') -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/LDF/ldfdyn.F90
r9598 r9923 408 408 zcmsmag = (rn_csmc/rpi)**2 ! (C_smag/pi)^2 409 409 zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag ) ! lower limit stability factor scaling 410 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * r dt )! upper limit stability factor scaling410 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rn_Dt ) ! upper limit stability factor scaling 411 411 IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo ! provide |U|L^3/12 lower limit instead 412 412 ! ! of |U|L^3/16 in blp case -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/LDF/ldftra.F90
r9737 r9923 852 852 ! 853 853 ! 854 zztmp = 0.5_wp * r au0 * rcp854 zztmp = 0.5_wp * rho0_rcp 855 855 IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN 856 856 zw2d(:,:) = 0._wp -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/OBS/diaobs.F90
r9656 r9923 539 539 ENDIF 540 540 541 idaystp = NINT( rday / r dt )541 idaystp = NINT( rday / rn_Dt ) 542 542 543 543 !----------------------------------------------------------------------- … … 630 630 631 631 ENDIF 632 632 ! 633 633 END SUBROUTINE dia_obs 634 634 635 635 636 SUBROUTINE dia_obs_wri … … 651 652 !! ! 15-08 (M. Martin) Combined writing for prof and surf types 652 653 !!---------------------------------------------------------------------- 653 !! * Modules used654 654 USE obs_rot_vel ! Rotation of velocities 655 655 656 656 IMPLICIT NONE 657 657 658 !! * Local declarations659 658 INTEGER :: jtype ! Data set loop variable 660 659 INTEGER :: jo, jvar, jk 661 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 662 & zu, & 663 & zv 660 REAL(wp), DIMENSION(:), ALLOCATABLE :: zu, zv 664 661 665 662 !----------------------------------------------------------------------- … … 771 768 !! ! 2014-09 (D. Lea) New generic routine now deals with arbitrary initial time of day 772 769 !!---------------------------------------------------------------------- 773 USE phycst, ONLY : & ! Physical constants 774 & rday 775 USE dom_oce, ONLY : & ! Ocean space and time domain variables 776 & rdt 770 USE phycst , ONLY : rday ! Physical constants 771 USE dom_oce, ONLY : rn_Dt ! Ocean space and time domain variables 777 772 778 773 IMPLICIT NONE 779 774 780 !! * Arguments 781 REAL(KIND=dp), INTENT(OUT) :: ddobs ! Date in YYYYMMDD.HHMMSS 782 INTEGER :: kstp 783 784 !! * Local declarations 775 REAL(KIND=dp), INTENT( out) :: ddobs ! Date in YYYYMMDD.HHMMSS 776 INTEGER , INTENT(in ) :: kstp 777 785 778 INTEGER :: iyea ! date - (year, month, day, hour, minute) 786 779 INTEGER :: imon … … 805 798 !! Compute number of days + number of hours + min since initial time 806 799 !!---------------------------------------------------------------------- 807 zdayfrc = kstp * r dt / rday800 zdayfrc = kstp * rn_Dt / rday 808 801 zdayfrc = zdayfrc - aint(zdayfrc) 809 802 imin = imin + int( zdayfrc * 24 * 60 ) … … 816 809 iday=iday+1 817 810 END DO 818 iday = iday + kstp * r dt / rday811 iday = iday + kstp * rn_Dt / rday 819 812 820 813 !----------------------------------------------------------------------- … … 842 835 END SUBROUTINE calc_date 843 836 837 844 838 SUBROUTINE ini_date( ddobsini ) 845 839 !!---------------------------------------------------------------------- … … 859 853 !! ! 2014-09 (D. Lea) Change to call generic routine calc_date 860 854 !!---------------------------------------------------------------------- 861 862 855 IMPLICIT NONE 863 864 !! * Arguments865 REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS866 856 ! 857 REAL(KIND=dp), INTENT(out) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS 858 !!---------------------------------------------------------------------- 859 ! 867 860 CALL calc_date( nit000 - 1, ddobsini ) 868 861 ! 869 862 END SUBROUTINE ini_date 863 870 864 871 865 SUBROUTINE fin_date( ddobsfin ) … … 1011 1005 END SUBROUTINE obs_setinterpopts 1012 1006 1007 !!====================================================================== 1013 1008 END MODULE diaobs -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/OBS/obs_prep.F90
r9598 r9923 610 610 !! ! 2010-05 (D. Lea) Fix in leap year calculation for NEMO vn3.2 611 611 !!---------------------------------------------------------------------- 612 !! * Modules used 613 USE dom_oce, ONLY : & ! Geographical information 614 & rdt 615 USE phycst, ONLY : & ! Physical constants 616 & rday, & 617 & rmmss, & 618 & rhhmm 619 !! * Arguments 612 USE dom_oce, ONLY : rn_Dt ! Geographical information 613 USE phycst , ONLY : rday, rmmss, rhhmm ! Physical constants 614 620 615 INTEGER, INTENT(IN) :: kcycle ! Current cycle 621 616 INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates … … 632 627 & kobshou, & 633 628 & kobsmin 634 INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 635 & kobsqc ! Quality control flag 636 INTEGER, DIMENSION(kobsno), INTENT(OUT) :: & 637 & kobsstp ! Number of time steps up to the 638 ! observation time 639 640 !! * Local declarations 629 INTEGER, DIMENSION(kobsno), INTENT(inout) :: kobsqc ! Quality control flag 630 INTEGER, DIMENSION(kobsno), INTENT( out) :: kobsstp ! Number of time steps up to the observation time 631 ! 641 632 INTEGER :: jyea 642 633 INTEGER :: jmon … … 661 652 662 653 ! Intialize the number of time steps per day 663 idaystp = NINT( rday / r dt )654 idaystp = NINT( rday / rn_Dt ) 664 655 665 656 !--------------------------------------------------------------------- … … 731 722 732 723 ! Add in the number of time steps to the observation minute 733 zminstp = rmmss / r dt724 zminstp = rmmss / rn_Dt 734 725 zhoustp = rhhmm * zminstp 735 726 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/fldread.F90
r9807 r9923 180 180 ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 181 181 IF( present(kit) ) THEN ! ignore kn_fsbc in this case 182 isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( r dt/REAL(nn_baro,wp) )182 isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rn_Dt/REAL(nn_e,wp) ) 183 183 ELSE ! middle of sbc time step 184 isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * r dt) + it_offset * NINT(rdt)184 isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rn_Dt) + it_offset * NINT(rn_Dt) 185 185 ENDIF 186 186 imf = SIZE( sd ) … … 213 213 CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit ) ! update after record informations 214 214 215 ! if kn_fsbc*r dt is larger than nfreqh (which is kind of odd),215 ! if kn_fsbc*rn_Dt is larger than nfreqh (which is kind of odd), 216 216 ! it is possible that the before value is no more the good one... we have to re-read it 217 217 ! if before is not the last record of the file currently opened and after is the first record to be read … … 234 234 IF( sd(jf)%ln_tint ) THEN 235 235 236 ! if kn_fsbc*r dt is larger than nfreqh (which is kind of odd),236 ! if kn_fsbc*rn_Dt is larger than nfreqh (which is kind of odd), 237 237 ! it is possible that the before value is no more the good one... we have to re-read it 238 238 ! if before record is not just just before the after record... … … 267 267 ! year/month/week/day file to be not present. If the run continue further than the current 268 268 ! year/month/week/day, next year/month/week/day file must exist 269 isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(r dt) ! second at the end of the run270 llstop = isecend > sd(jf)%nrec_a(2) 269 isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rn_Dt) ! second at the end of the run 270 llstop = isecend > sd(jf)%nrec_a(2) ! read more than 1 record of next year 271 271 ! we suppose that the date of next file is next day (should be ok even for weekly files...) 272 272 CALL fld_clopn( sd(jf), nyear + COUNT((/llnxtyr /)) , & … … 485 485 ENDIF 486 486 IF( PRESENT(kt_offset) ) it_offset = kt_offset 487 IF( PRESENT(kit) ) THEN ; it_offset = ( kit + it_offset ) * NINT( r dt/REAL(nn_baro,wp) )488 ELSE ; it_offset = it_offset * NINT( rdt)487 IF( PRESENT(kit) ) THEN ; it_offset = ( kit + it_offset ) * NINT( rn_Dt / REAL(nn_e,wp) ) 488 ELSE ; it_offset = it_offset * NINT( rn_Dt ) 489 489 ENDIF 490 490 ! … … 563 563 ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year 564 564 ENDIF 565 ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * r dt + REAL( it_offset, wp )! centrered in the middle of sbc time step566 ztmp = ztmp + 0.01 * r dt! avoid truncation error565 ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rn_Dt + REAL( it_offset, wp ) ! centrered in the middle of sbc time step 566 ztmp = ztmp + 0.01 * rn_Dt ! avoid truncation error 567 567 IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record 568 568 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcapr.F90
r9598 r9923 36 36 37 37 REAL(wp) :: tarea ! whole domain mean masked ocean surface 38 REAL(wp) :: r1_ grau ! = 1.e0 / (grav * rau0)38 REAL(wp) :: r1_rhog ! = 1 / (rho0*grav) 39 39 40 40 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_apr ! structure of input fields (file informations, fields read) … … 100 100 ENDIF 101 101 ! 102 r1_ grau = 1.e0 / (grav * rau0)!* constant for optimization102 r1_rhog = 1._wp / (rho0*grav) !* constant for optimization 103 103 ! 104 104 ! !* control check … … 144 144 ! 145 145 ! !* Patm related forcing at kt 146 ssh_ib(:,:) = - ( sf_apr(1)%fnow(:,:,1) - rn_pref ) * r1_ grau! equivalent ssh (inverse barometer)146 ssh_ib(:,:) = - ( sf_apr(1)%fnow(:,:,1) - rn_pref ) * r1_rhog ! equivalent ssh (inverse barometer) 147 147 apr (:,:) = sf_apr(1)%fnow(:,:,1) ! atmospheric pressure 148 148 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcblk.F90
r9767 r9923 225 225 ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 226 226 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 227 IF( slf_i(ifpr)%nfreqh > 0. .AND. MOD( 3600. * slf_i(ifpr)%nfreqh , REAL(nn_fsbc) * r dt) /= 0. ) &228 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep r dt*nn_fsbc is NOT a submultiple of atmosphericforcing frequency.', &229 & ' This is not ideal. You should consider changing either r dt or nn_fsbc value...' )227 IF( slf_i(ifpr)%nfreqh > 0. .AND. MOD( 3600. * slf_i(ifpr)%nfreqh , REAL(nn_fsbc) * rn_Dt) /= 0. ) & 228 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmos. forcing frequency.', & 229 & ' This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' ) 230 230 231 231 END DO … … 323 323 ! 324 324 ! ! compute the surface ocean fluxes using bulk formulea 325 IF( MOD( kt -1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )325 IF( MOD( kt-1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 326 326 327 327 #if defined key_cice 328 IF( MOD( kt -1, nn_fsbc ) == 0 ) THEN328 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 329 329 qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) 330 330 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbccpl.F90
r9767 r9923 193 193 194 194 REAL(wp) :: rpref = 101000._wp ! reference atmospheric pressure[N/m2] 195 REAL(wp) :: r1_ grau ! = 1.e0 / (grav * rau0)195 REAL(wp) :: r1_rhog ! = 1 / (rho0*grav) 196 196 197 197 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) :: nrcvinfo ! OASIS info argument … … 1100 1100 LOGICAL :: llnewtx, llnewtau ! update wind stress components and module?? 1101 1101 INTEGER :: ji, jj, jn ! dummy loop indices 1102 INTEGER :: isec ! number of seconds since nit000 (assuming r dt did not change since nit000)1102 INTEGER :: isec ! number of seconds since nit000 (assuming rn_Dt did not change since nit000) 1103 1103 REAL(wp) :: zcumulneg, zcumulpos ! temporary scalars 1104 1104 REAL(wp) :: zcoef ! temporary scalar … … 1114 1114 ! ! Receive all the atmos. fields (including ice information) 1115 1115 ! ! ======================================================= ! 1116 isec = ( kt - nit000 ) * NINT( r dt )! date of exchanges1116 isec = ( kt - nit000 ) * NINT( rn_Dt ) ! date of exchanges 1117 1117 DO jn = 1, jprcv ! received fields sent by the atmosphere 1118 1118 IF( srcv(jn)%laction ) CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) ) … … 1259 1259 IF( kt /= nit000 ) ssh_ibb(:,:) = ssh_ib(:,:) !* Swap of ssh_ib fields 1260 1260 1261 r1_ grau = 1.e0 / (grav * rau0) !* constant for optimization1262 ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_ grau! equivalent ssh (inverse barometer)1261 r1_rhog = 1.e0 / (grav * rho0) !* constant for optimization 1262 ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_rhog ! equivalent ssh (inverse barometer) 1263 1263 apr (:,:) = frcv(jpr_mslp)%z3(:,:,1) !atmospheric pressure 1264 1264 … … 2047 2047 !!---------------------------------------------------------------------- 2048 2048 ! 2049 isec = ( kt - nit000 ) * NINT( r dt ) ! date of exchanges2049 isec = ( kt - nit000 ) * NINT( rn_Dt ) ! date of exchanges 2050 2050 2051 2051 zfr_l(:,:) = 1.- fr_i(:,:) -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcdcy.F90
r9598 r9923 88 88 89 89 ! When are we during the day (from 0 to 1) 90 zlo = ( REAL(nsec_day, wp) - 0.5_wp * r dt ) / rday91 zup = zlo + ( REAL(nn_fsbc, wp) * r dt ) / rday90 zlo = ( REAL(nsec_day, wp) - 0.5_wp * rn_Dt ) / rday 91 zup = zlo + ( REAL(nn_fsbc, wp) * rn_Dt ) / rday 92 92 ! 93 93 IF( nday_qsr == -1 ) THEN ! first time step only … … 187 187 END DO 188 188 ! 189 ztmp = rday / ( r dt * REAL(nn_fsbc, wp) )189 ztmp = rday / ( rn_Dt * REAL(nn_fsbc, wp) ) 190 190 rscal(:,:) = rscal(:,:) * ztmp 191 191 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcfwb.F90
r9598 r9923 123 123 ENDIF 124 124 ! ! Update fwfold if new year start 125 ikty = 365 * 86400 / r dt!!bug use of 365 days leap year or 360d year !!!!!!!125 ikty = 365 * 86400 / rn_Dt !!bug use of 365 days leap year or 360d year !!!!!!! 126 126 IF( MOD( kt, ikty ) == 0 ) THEN 127 127 a_fwb_b = a_fwb ! mean sea level taking into account the ice+snow 128 128 ! sum over the global domain 129 a_fwb = glob_sum( e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_r au0 ) )129 a_fwb = glob_sum( e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rho0 ) ) 130 130 a_fwb = a_fwb * 1.e+3 / ( area * rday * 365. ) ! convert in Kg/m3/s = mm/s 131 131 !!gm ! !!bug 365d year -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcice_cice.F90
r9598 r9923 13 13 USE dom_oce ! ocean space and time domain 14 14 USE domvvl 15 USE phycst , only : rcp, rau0, r1_rau0, rhosn, rhoic15 USE phycst , ONLY : rcp, rho0, r1_rho0, rhosn, rhoic 16 16 USE in_out_manager ! I/O manager 17 17 USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit … … 227 227 IF( .NOT.ln_rstart ) THEN 228 228 IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area 229 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_r au0230 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_r au0229 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rho0 230 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rho0 231 231 232 232 !!gm This should be put elsewhere.... (same remark for limsbc) … … 422 422 ! Freezing/melting potential 423 423 ! Calculated over NEMO leapfrog timestep (hence 2*dt) 424 nfrzmlt(:,:) = r au0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt )424 nfrzmlt(:,:) = rho0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt ) 425 425 426 426 ztmp(:,:) = nfrzmlt(:,:) … … 459 459 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 460 460 ! 461 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_r au0461 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rho0 462 462 ! 463 463 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcisf.F90
r9728 r9923 54 54 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp !: specific heat of ice shelf [J/kg/K] 55 55 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp !: heat diffusivity through the ice-shelf [m2/s] 56 REAL(wp), PUBLIC, SAVE :: rho isf= 920.0_wp !: volumic mass of ice shelf [kg/m3]56 REAL(wp), PUBLIC, SAVE :: rho_isf = 920.0_wp !: volumic mass of ice shelf [kg/m3] 57 57 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp !: air temperature on top of ice shelf [C] 58 58 REAL(wp), PUBLIC, SAVE :: rlfusisf = 0.334e6_wp !: latent heat of fusion of ice shelf [J/kg] … … 144 144 ! compute tsc due to isf 145 145 ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 146 ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / r au0).146 ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rho0). 147 147 ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 148 148 DO jj = 1,jpj … … 153 153 CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 154 154 155 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_r au0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !155 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rho0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rho0 ! 156 156 risf_tsc(:,:,jp_sal) = 0.0_wp 157 157 … … 160 160 ! output 161 161 IF( iom_use('iceshelf_cea') ) CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) ) ! isf mass flux 162 IF( iom_use('hflx_isf_cea') ) CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * r au0 * rcp ) ! isf sensible+latent heat (W/m2)162 IF( iom_use('hflx_isf_cea') ) CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rho0 * rcp ) ! isf sensible+latent heat (W/m2) 163 163 IF( iom_use('qlatisf' ) ) CALL iom_put( 'qlatisf' , qisf(:,:) ) ! isf latent heat 164 164 IF( iom_use('fwfisf' ) ) CALL iom_put( 'fwfisf' , fwfisf(:,:) ) ! isf mass flux (opposite sign) … … 452 452 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 453 453 ! For those corresponding to zonal boundary 454 qisf(ji,jj) = - r au0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave &454 qisf(ji,jj) = - rho0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 455 455 & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 456 456 … … 500 500 zlamb1 =-0.0564_wp 501 501 zlamb2 = 0.0773_wp 502 zlamb3 =-7.8633e-8 * grav * r au0502 zlamb3 =-7.8633e-8 * grav * rho0 503 503 ELSE ! linearisation from table 4 (Asay-Davis et al., 2015) 504 504 zlamb1 =-0.0573_wp 505 505 zlamb2 = 0.0832_wp 506 zlamb3 =-7.53e-8 * grav * r au0506 zlamb3 =-7.53e-8 * grav * rho0 507 507 ENDIF 508 508 ! … … 526 526 DO jj = 1, jpj 527 527 DO ji = 1, jpi 528 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*r au0*(ttbl(ji,jj)-zfrz(ji,jj))528 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rho0*(ttbl(ji,jj)-zfrz(ji,jj)) 529 529 zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 530 530 END DO … … 544 544 DO ji = 1, jpi 545 545 ! compute coeficient to solve the 2nd order equation 546 zeps1 = rcp*r au0*zgammat(ji,jj)547 zeps2 = rlfusisf*r au0*zgammas(ji,jj)548 zeps3 = rho isf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)546 zeps1 = rcp*rho0*zgammat(ji,jj) 547 zeps2 = rlfusisf*rho0*zgammas(ji,jj) 548 zeps3 = rho_isf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 549 549 zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 550 550 zeps6 = zeps4-ttbl(ji,jj) … … 567 567 ! zhtflx is upward heat flux (out of ocean) 568 568 ! compute the upward water and heat flux (eq. 28 and eq. 29) 569 zfwflx(ji,jj) = r au0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)570 zhtflx(ji,jj) = zgammat(ji,jj) * r au0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) )569 zfwflx(ji,jj) = rho0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 570 zhtflx(ji,jj) = zgammat(ji,jj) * rho0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 571 571 END DO 572 572 END DO … … 890 890 DO jk = ikt, ikb - 1 891 891 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 892 & * r1_hisf_tbl(ji,jj) * r1_r au0 * zfact892 & * r1_hisf_tbl(ji,jj) * r1_rho0 * zfact 893 893 END DO 894 894 ! level partially include in ice shelf boundary layer 895 895 phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 896 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_r au0 * zfact * ralpha(ji,jj)896 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rho0 * zfact * ralpha(ji,jj) 897 897 END DO 898 898 END DO -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcmod.F90
r9656 r9923 177 177 ! 178 178 IF( .NOT.ln_usr ) THEN ! the model calendar needs some specificities (except in user defined case) 179 IF( MOD( rday , rdt ) /= 0. ) CALL ctl_stop( 'the time step must devide the number of second of in a day' )180 IF( MOD( rday ,2. ) /= 0. ) CALL ctl_stop( 'the number of second of in a day must be an even number' )181 IF( MOD( r dt , 2.) /= 0. ) CALL ctl_stop( 'the time step (in second) must be an even number' )179 IF( MOD( rday , rn_Dt ) /= 0. ) CALL ctl_stop( 'the time step must devide the number of second of in a day' ) 180 IF( MOD( rday , 2. ) /= 0. ) CALL ctl_stop( 'the number of second of in a day must be an even number' ) 181 IF( MOD( rn_Dt, 2. ) /= 0. ) CALL ctl_stop( 'the time step (in second) must be an even number' ) 182 182 ENDIF 183 183 ! !** check option consistency … … 288 288 ! SAS time-step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly 289 289 IF( nn_components /= jp_iam_nemo ) THEN 290 IF( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(r dt)291 IF( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(r dt)290 IF( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rn_Dt) 291 IF( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rn_Dt) 292 292 ! 293 293 IF(lwp)THEN … … 306 306 ENDIF 307 307 ! 308 IF( MOD( rday, REAL(nn_fsbc, wp) * r dt ) /= 0 ) &308 IF( MOD( rday, REAL(nn_fsbc, wp) * rn_Dt ) /= 0 ) & 309 309 & CALL ctl_warn( 'sbc_init : nn_fsbc is NOT a multiple of the number of time steps in a day' ) 310 310 ! 311 IF( ln_dm2dc .AND. NINT(rday) / ( nn_fsbc * NINT(r dt) ) < 8 ) &311 IF( ln_dm2dc .AND. NINT(rday) / ( nn_fsbc * NINT(rn_Dt) ) < 8 ) & 312 312 & CALL ctl_warn( 'sbc_init : diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' ) 313 313 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcrnf.F90
r9727 r9923 116 116 IF( ln_rnf_sal ) CALL fld_read ( kt, nn_fsbc, sf_s_rnf ) ! idem for runoffs salinity if required 117 117 ! 118 IF( MOD( kt -1, nn_fsbc ) == 0 ) THEN118 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 119 119 ! 120 120 IF( .NOT. l_rnfcpl ) rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1) ! updated runoff value at time step kt … … 122 122 ! ! set temperature & salinity content of runoffs 123 123 IF( ln_rnf_tem ) THEN ! use runoffs temperature data 124 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_r au0124 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0 125 125 CALL eos_fzp( sss_m(:,:), ztfrz(:,:) ) 126 126 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp ) ! if missing data value use SST as runoffs temperature 127 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_r au0127 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rho0 128 128 END WHERE 129 129 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp ) ! where fwf comes from melting of ice shelves or iceberg 130 rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_r au0 - rnf(:,:) * rlfusisf * r1_rau0_rcp130 rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rho0 - rnf(:,:) * rlfusisf * r1_rho0_rcp 131 131 END WHERE 132 132 ELSE ! use SST as runoffs temperature 133 133 !CEOD River is fresh water so must at least be 0 unless we consider ice 134 rnf_tsc(:,:,jp_tem) = MAX(sst_m(:,:),0.0_wp) * rnf(:,:) * r1_r au0134 rnf_tsc(:,:,jp_tem) = MAX(sst_m(:,:),0.0_wp) * rnf(:,:) * r1_rho0 135 135 ENDIF 136 136 ! ! use runoffs salinity data 137 IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_r au0137 IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0 138 138 ! ! else use S=0 for runoffs (done one for all in the init) 139 139 IF( iom_use('runoffs') ) CALL iom_put( 'runoffs' , rnf(:,:) ) ! output runoff mass flux 140 IF( iom_use('hflx_rnf_cea') ) CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * r au0 * rcp ) ! output runoff sensible heat (W/m2)140 IF( iom_use('hflx_rnf_cea') ) CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rho0 * rcp ) ! output runoff sensible heat (W/m2) 141 141 ENDIF 142 142 ! … … 198 198 DO ji = 1, jpi 199 199 DO jk = 1, nk_rnf(ji,jj) 200 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_r au0 / h_rnf(ji,jj)200 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj) 201 201 END DO 202 202 END DO … … 211 211 ! ! apply the runoff input flow 212 212 DO jk = 1, nk_rnf(ji,jj) 213 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_r au0 / h_rnf(ji,jj)213 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj) 214 214 END DO 215 215 END DO … … 218 218 ELSE !== runoff put only at the surface ==! 219 219 h_rnf (:,:) = e3t_n (:,:,1) ! update h_rnf to be depth of top box 220 phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_r au0 / e3t_n(:,:,1)220 phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rho0 / e3t_n(:,:,1) 221 221 ENDIF 222 222 ! -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcssm.F90
r9598 r9923 106 106 frq_m(:,:) = zcoef * fraqsr_1lev(:,:) 107 107 ! ! ---------------------------------------- ! 108 ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN! Initialisation: New mean computation !108 ELSEIF( MOD( kt-2, nn_fsbc ) == 0 ) THEN ! Initialisation: New mean computation ! 109 109 ! ! ---------------------------------------- ! 110 110 ssu_m(:,:) = 0._wp ! reset to zero ocean mean sbc fields … … 135 135 136 136 ! ! ---------------------------------------- ! 137 IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN! Mean value at each nn_fsbc time-step !137 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Mean value at each nn_fsbc time-step ! 138 138 ! ! ---------------------------------------- ! 139 139 zcoef = 1. / REAL( nn_fsbc, wp ) … … 263 263 CALL iom_set_rstw_var_active('frq_m') 264 264 ENDIF 265 265 ! 266 266 END SUBROUTINE sbc_ssm_init 267 267 -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbctide.F90
r9598 r9923 48 48 !!---------------------------------------------------------------------- 49 49 50 IF( nsec_day == NINT(0.5_wp * r dt) .OR. kt == nit000 ) THEN ! start a new day50 IF( nsec_day == NINT(0.5_wp * rn_Dt) .OR. kt == nit000 ) THEN ! start a new day 51 51 ! 52 52 IF( kt == nit000 )THEN … … 72 72 ! Temporarily set nsec_day to beginning of day. 73 73 nsec_day_orig = nsec_day 74 IF ( nsec_day /= NINT(0.5_wp * r dt) ) THEN75 kt_tide = kt - (nsec_day - 0.5_wp * r dt)/rdt76 nsec_day = NINT(0.5_wp * r dt)74 IF ( nsec_day /= NINT(0.5_wp * rn_Dt) ) THEN 75 kt_tide = kt - (nsec_day - 0.5_wp * rn_Dt) / rn_Dt 76 nsec_day = NINT(0.5_wp * rn_Dt) 77 77 ELSE 78 78 kt_tide = kt -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/tideini.F90
r9598 r9923 20 20 PUBLIC 21 21 22 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: omega_tide !: 23 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: v0tide !: 24 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: utide !: 25 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: ftide !: 22 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: omega_tide !: 23 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: v0tide !: 24 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: utide !: 25 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: ftide !: 26 26 27 LOGICAL , PUBLIC :: ln_tide !: 28 LOGICAL , PUBLIC :: ln_tide_pot !: 29 LOGICAL , PUBLIC :: ln_read_load !: 30 LOGICAL , PUBLIC :: ln_scal_load !: 31 LOGICAL , PUBLIC :: ln_tide_ramp !: 32 INTEGER , PUBLIC :: nb_harmo !: 33 INTEGER , PUBLIC :: kt_tide !: 34 REAL(wp), PUBLIC :: rdttideramp !: 35 REAL(wp), PUBLIC :: rn_scal_load !: 36 CHARACTER(lc), PUBLIC :: cn_tide_load !: 27 ! !!* nam_tide namelist * 28 LOGICAL , PUBLIC :: ln_tide !: Use tidal components 29 LOGICAL , PUBLIC :: ln_tide_pot !: Apply astronomical potential 30 LOGICAL , PUBLIC :: ln_read_load !: Read load potential from file 31 CHARACTER(lc), PUBLIC :: cn_tide_load !: associated file name 32 LOGICAL , PUBLIC :: ln_scal_load !: Use a scalar approximation for load potential 33 REAL(wp), PUBLIC :: rn_load !: SSH fraction used in scalar approximation 34 LOGICAL , PUBLIC :: ln_tide_ramp !: Apply ramp on tides at startup 35 REAL(wp), PUBLIC :: rn_ramp !: Duration of ramp [days] 36 INTEGER , PUBLIC :: nb_harmo !: number of tidal harmonique used 37 INTEGER , PUBLIC :: kt_tide !: ??? 37 38 38 INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide !: 39 INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide !: ??? 39 40 40 41 !!---------------------------------------------------------------------- … … 52 53 CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 53 54 INTEGER :: ios ! Local integer output status for namelist read 54 ! 55 !! 55 56 NAMELIST/nam_tide/ln_tide, ln_tide_pot, ln_scal_load, ln_read_load, cn_tide_load, & 56 & ln_tide_ramp, rn_ scal_load, rdttideramp, clname57 & ln_tide_ramp, rn_load, rn_ramp, clname 57 58 !!---------------------------------------------------------------------- 58 59 ! … … 76 77 WRITE(numout,*) ' Apply astronomical potential ln_tide_pot = ', ln_tide_pot 77 78 WRITE(numout,*) ' Use scalar approx. for load potential ln_scal_load = ', ln_scal_load 79 WRITE(numout,*) ' SSH fraction used in scal. approx. rn_load = ', rn_load 78 80 WRITE(numout,*) ' Read load potential from file ln_read_load = ', ln_read_load 79 81 WRITE(numout,*) ' Apply ramp on tides at startup ln_tide_ramp = ', ln_tide_ramp 80 WRITE(numout,*) ' Fraction of SSH used in scal. approx. rn_scal_load = ', rn_scal_load 81 WRITE(numout,*) ' Duration (days) of ramp rdttideramp = ', rdttideramp 82 WRITE(numout,*) ' Duration of ramp rn_ramp = ', rn_ramp, ' [days]' 82 83 ENDIF 83 84 ELSE 84 rn_ scal_load = 0._wp85 85 rn_load = 0._wp 86 ! 86 87 IF(lwp) WRITE(numout,*) 87 88 IF(lwp) WRITE(numout,*) 'tide_init : tidal components not used (ln_tide = F)' … … 92 93 CALL tide_init_Wave 93 94 ! 94 nb_harmo =095 nb_harmo = 0 95 96 DO jk = 1, jpmax_harmo 96 97 DO ji = 1,jpmax_harmo … … 108 109 IF( ln_scal_load.AND.ln_read_load ) & 109 110 & CALL ctl_stop('Choose between ln_scal_load and ln_read_load') 110 IF( ln_tide_ramp.AND.((nitend-nit000+1)*r dt/rday < rdttideramp) ) &111 & CALL ctl_stop('r dttideramp must be lower than run duration')112 IF( ln_tide_ramp.AND.(r dttideramp<0.) ) &113 & CALL ctl_stop('r dttideramp must be positive')111 IF( ln_tide_ramp.AND.((nitend-nit000+1)*rn_Dt/rday < rn_ramp) ) & 112 & CALL ctl_stop('rn_ramp must be lower than run duration') 113 IF( ln_tide_ramp.AND.(rn_ramp<0.) ) & 114 & CALL ctl_stop('rn_ramp must be positive') 114 115 ! 115 116 ALLOCATE( ntide(nb_harmo) ) … … 123 124 END DO 124 125 ! 125 ALLOCATE( omega_tide(nb_harmo), v0tide 126 & utide (nb_harmo), ftide 126 ALLOCATE( omega_tide(nb_harmo), v0tide(nb_harmo), & 127 & utide (nb_harmo), ftide (nb_harmo) ) 127 128 kt_tide = nit000 128 129 ! 129 IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp130 IF (.NOT.ln_scal_load ) rn_load = 0._wp 130 131 ! 131 132 END SUBROUTINE tide_init -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/updtide.F90
r9598 r9923 6 6 !! History : 9.0 ! 07 (O. Le Galloudec) Original code 7 7 !!---------------------------------------------------------------------- 8 !! upd_tide : update tidal potential 8 9 9 !!---------------------------------------------------------------------- 10 USE oce ! ocean dynamics and tracers variables 11 USE dom_oce ! ocean space and time domain 12 USE in_out_manager ! I/O units 13 USE phycst ! physical constant 14 USE sbctide ! tide potential variable 15 USE tideini, ONLY: ln_tide_ramp, rdttideramp 10 !! upd_tide : update tidal potential 11 !!---------------------------------------------------------------------- 12 USE oce ! ocean dynamics and tracers variables 13 USE dom_oce ! ocean space and time domain 14 USE in_out_manager ! I/O units 15 USE phycst ! physical constant 16 USE sbctide ! tide potential variable 17 USE tideini , ONLY : ln_tide_ramp, rn_ramp 16 18 17 19 IMPLICIT NONE … … 37 39 !! ** Action : pot_astro actronomical potential 38 40 !!---------------------------------------------------------------------- 39 INTEGER, INTENT(in) :: kt ! ocean time-step index 40 INTEGER, INTENT(in), OPTIONAL :: kit ! external mode sub-time-step index (lk_dynspg_ts=T) 41 INTEGER, INTENT(in), OPTIONAL :: time_offset ! time offset in number 42 ! of internal steps (lk_dynspg_ts=F) 43 ! of external steps (lk_dynspg_ts=T) 41 INTEGER, INTENT(in) :: kt ! ocean time-step index 42 INTEGER, INTENT(in), OPTIONAL :: kit ! external mode sub-time-step index (lk_dynspg_ts=T) 43 INTEGER, INTENT(in), OPTIONAL :: time_offset ! time offset in number of internal steps (lk_dynspg_ts=F) 44 ! ! of external steps (lk_dynspg_ts=T) 44 45 ! 46 INTEGER :: ji, jj, jk ! dummy loop indices 45 47 INTEGER :: joffset ! local integer 46 INTEGER :: ji, jj, jk ! dummy loop indices47 48 REAL(wp) :: zt, zramp ! local scalar 48 49 REAL(wp), DIMENSION(nb_harmo) :: zwt … … 50 51 ! 51 52 ! ! tide pulsation at model time step (or sub-time-step) 52 zt = ( kt - kt_tide ) * r dt53 zt = ( kt - kt_tide ) * rn_Dt 53 54 ! 54 55 joffset = 0 … … 56 57 ! 57 58 IF( PRESENT( kit ) ) THEN 58 zt = zt + ( kit + joffset - 1 ) * r dt / REAL( nn_baro, wp )59 zt = zt + ( kit + joffset - 1 ) * rn_Dt / REAL( nn_e, wp ) 59 60 ELSE 60 zt = zt + joffset * r dt61 zt = zt + joffset * rn_Dt 61 62 ENDIF 62 63 ! … … 69 70 ! 70 71 IF( ln_tide_ramp ) THEN ! linear increase if asked 71 zt = ( kt - nit000 ) * r dt72 IF( PRESENT( kit ) ) zt = zt + ( kit + joffset -1) * r dt / REAL( nn_baro, wp )73 zramp = MIN( MAX( zt / (rdttideramp*rday) , 0._wp) , 1._wp )72 zt = ( kt - nit000 ) * rn_Dt 73 IF( PRESENT( kit ) ) zt = zt + ( kit + joffset -1) * rn_Dt / REAL( nn_e, wp ) 74 zramp = MIN( MAX( 0._wp , zt / (rn_ramp*rday) ) , 1._wp ) 74 75 pot_astro(:,:) = zramp * pot_astro(:,:) 75 76 ENDIF -
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/eosbn2.F90