Changeset 13877
- Timestamp:
- 2020-11-25T15:23:00+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3
- Files:
-
- 70 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette@13 292sette10 ^/utils/CI/sette@13559 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/INSTALL.rst
r11734 r13877 122 122 .. code:: console 123 123 124 $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/ trunk124 $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/r4.0/r4.0.4 125 125 126 126 Description of 1\ :sup:`st` level tree structure -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-ice.xml
r9896 r13877 52 52 <field field_ref="normstr" name="normstr" /> 53 53 <field field_ref="sheastr" name="sheastr" /> 54 <field field_ref="isig1" name="isig1" /> 55 <field field_ref="isig2" name="isig2" /> 56 <field field_ref="isig3" name="isig3" /> 54 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 55 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 57 56 58 57 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/C1D_PAPA/EXPREF/namelist_cfg
r13631 r13877 428 428 / 429 429 !----------------------------------------------------------------------- 430 !-----------------------------------------------------------------------431 /432 !-----------------------------------------------------------------------433 430 &namhsb ! Heat and salt budgets (default: OFF) 434 431 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-ice.xml
r11945 r13877 53 53 <field field_ref="normstr" name="normstr" /> 54 54 <field field_ref="sheastr" name="sheastr" /> 55 <field field_ref="isig1" name="isig1" /> 56 <field field_ref="isig2" name="isig2" /> 57 <field field_ref="isig3" name="isig3" /> 55 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 56 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 58 57 59 58 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml
r12377 r13877 53 53 <field field_ref="normstr" name="normstr" /> 54 54 <field field_ref="sheastr" name="sheastr" /> 55 <field field_ref="isig1" name="isig1" /> 56 <field field_ref="isig2" name="isig2" /> 57 <field field_ref="isig3" name="isig3" /> 55 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 56 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 58 57 59 58 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg
r13631 r13877 388 388 / 389 389 !----------------------------------------------------------------------- 390 !-----------------------------------------------------------------------391 /392 !-----------------------------------------------------------------------393 390 &namhsb ! Heat and salt budgets (default: OFF) 394 391 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg
r13631 r13877 384 384 / 385 385 !----------------------------------------------------------------------- 386 !-----------------------------------------------------------------------387 /388 !-----------------------------------------------------------------------389 386 &namhsb ! Heat and salt budgets (default: OFF) 390 387 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml
r9896 r13877 52 52 <field field_ref="normstr" name="normstr" /> 53 53 <field field_ref="sheastr" name="sheastr" /> 54 <field field_ref="isig1" name="isig1" /> 55 <field field_ref="isig2" name="isig2" /> 56 <field field_ref="isig3" name="isig3" /> 54 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 55 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 57 56 58 57 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/field_def_nemo-ice.xml
r13631 r13877 75 75 <field id="utau_bi" long_name="X-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_x_stress" unit="N/m2" /> 76 76 <field id="vtau_bi" long_name="Y-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_y_stress" unit="N/m2" /> 77 <field id="isig1" long_name="1st principal stress component for EVP rhg" unit="" /> 78 <field id="isig2" long_name="2nd principal stress component for EVP rhg" unit="" /> 79 <field id="isig3" long_name="convergence measure for EVP rheology (must be around 1)" unit="" /> 77 <field id="sig1_pnorm" long_name="P-normalized 1st principal stress component" unit="" /> 78 <field id="sig2_pnorm" long_name="P-normalized 2nd principal stress component" unit="" /> 80 79 <field id="normstr" long_name="Average normal stress in sea ice" standard_name="average_normal_stress" unit="N/m" /> 81 80 <field id="sheastr" long_name="Maximum shear stress in sea ice" standard_name="maximum_shear_stress" unit="N/m" /> … … 165 164 166 165 <!-- diags --> 167 <field id="hfxdhc" long_name="Heat content variation in snow and ice (neg = ice cooling)" unit="W/m2" /> 166 <field id="hfxdhc" long_name="Heat content variation in snow and ice (neg = ice cooling)" unit="W/m2" /> 167 <!-- available if ln_icediachk=T --> 168 <field id="icedrift_mass" long_name="Ice mass drift (conservation check)" unit="kg/m2/s" /> 169 <field id="icedrift_salt" long_name="Ice salt drift (conservation check)" unit="kg/m2/s" /> 170 <field id="icedrift_heat" long_name="Ice heat drift (conservation check)" unit="W/m2" /> 168 171 169 172 <!-- sbcssm variables --> … … 400 403 <field field_ref="normstr" name="normstr" /> 401 404 <field field_ref="sheastr" name="sheastr" /> 402 <field field_ref="isig1" name="isig1" /> 403 <field field_ref="isig2" name="isig2" /> 404 <field field_ref="isig3" name="isig3" /> 405 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 406 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 405 407 406 408 <!-- heat fluxes --> … … 459 461 <field field_ref="vfxthin" name="vfxthin" /> 460 462 461 <!-- diag error for negative ice volume after advection -->462 <field field_ref="iceneg_pres" name="sineg_pres" />463 <field field_ref="iceneg_volu" name="sineg_volu" />464 <field field_ref="iceneg_hfx" name="sineg_hfx" />465 463 </field_group> 466 464 -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/namelist_ref
r13631 r13877 1264 1264 !!gm ln_trdmld_instant = .false. ! flag to diagnose trends of instantantaneous or mean ML T/S 1265 1265 !!gm 1266 /1267 1266 !----------------------------------------------------------------------- 1268 1267 &namhsb ! Heat and salt budgets (default: OFF) -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml
r11536 r13877 53 53 <field field_ref="normstr" name="normstr" /> 54 54 <field field_ref="sheastr" name="sheastr" /> 55 <field field_ref="isig1" name="isig1" /> 56 <field field_ref="isig2" name="isig2" /> 57 <field field_ref="isig3" name="isig3" /> 55 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 56 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 58 57 59 58 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/file_def_nemo-ice.xml
r12905 r13877 50 50 <field field_ref="normstr" name="normstr" /> 51 51 <field field_ref="sheastr" name="sheastr" /> 52 <field field_ref="isig1" name="isig1" /> 53 <field field_ref="isig2" name="isig2" /> 54 <field field_ref="isig3" name="isig3" /> 52 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 53 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 55 54 56 55 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/namelist_cfg
r13631 r13877 583 583 / 584 584 !----------------------------------------------------------------------- 585 !-----------------------------------------------------------------------586 /587 !-----------------------------------------------------------------------588 585 &namhsb ! Heat and salt budgets (default: OFF) 589 586 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/NEMO_manual_state.txt
r10569 r13877 39 39 namdia: iiceprt jiceprt 40 40 nam_diaharm: nit000_han nitend_han nstep_han tname(1) tname(2) 41 namdrg: ln_ OFF41 namdrg: ln_drg_OFF 42 42 namdrg_bot: rn_Cd0 rn_Uc0 rn_Cdmax 43 43 namdrg_top: rn_Cd0 rn_Uc0 rn_Cdmax -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_TRA.tex
r11693 r13877 1229 1229 In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, 1230 1230 with $\rho_o$ a reference density. 1231 Called \textit{r au0} in the code,1231 Called \textit{rho0} in the code, 1232 1232 $\rho_o$ is set in \mdl{phycst} to a value of \texttt{1,026} $Kg/m^3$. 1233 1233 This is a sensible choice for the reference density used in a Boussinesq ocean climate model, -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_ZDF.tex
r11693 r13877 1160 1160 \] 1161 1161 When \np[=.true.]{ln_lin}{ln\_lin}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}. 1162 Setting \np[=.true.]{ln_ OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.1162 Setting \np[=.true.]{ln_drg_OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition. 1163 1163 1164 1164 These values are assigned in \mdl{zdfdrg}. -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nambdy_dta
r11703 r13877 29 29 bn_aip = 'NOT USED' , 24. , 'siapnd' , .true. , .false., 'daily' , '' , '' , '' 30 30 bn_hip = 'NOT USED' , 24. , 'sihpnd' , .true. , .false., 'daily' , '' , '' , '' 31 bn_hil = 'NOT USED' , 24. , 'sihlid' , .true. , .false., 'daily' , '' , '' , '' 31 32 ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds 32 33 rn_ice_tem = 270. ! arbitrary temperature of incoming sea ice … … 35 36 rn_ice_apnd = 0.2 ! -- pond fraction = a_ip/a_i -- 36 37 rn_ice_hpnd = 0.05 ! -- pond depth -- 38 rn_ice_hlid = 0.0 ! -- pond lid depth -- 37 39 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdia
r11703 r13877 8 8 ln_icediahsb = .false. ! output the heat, mass & salt budgets (T) or not (F) 9 9 ln_icectl = .false. ! ice points output for debug (T or F) 10 iiceprt = 10 !i-index for debug11 jiceprt = 10 !j-index for debug10 iiceprt = 10 ! i-index for debug 11 jiceprt = 10 ! j-index for debug 12 12 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg
r10075 r13877 2 2 &namdrg ! top/bottom drag coefficient (default: NO selection) 3 3 !----------------------------------------------------------------------- 4 ln_ OFF= .false. ! free-slip : Cd = 0 (F => fill namdrg_bot4 ln_drg_OFF = .false. ! free-slip : Cd = 0 (F => fill namdrg_bot 5 5 ln_lin = .false. ! linear drag: Cd = Cd0 Uc0 & namdrg_top) 6 6 ln_non_lin = .false. ! non-linear drag: Cd = Cd0 |U| … … 8 8 ! 9 9 ln_drgimp = .true. ! implicit top/bottom friction flag 10 ln_drgice_imp = .false. ! implicit ice-ocean drag 10 11 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_bot
r10075 r13877 1 1 !----------------------------------------------------------------------- 2 &namdrg_bot ! BOTTOM friction (ln_ OFF =F)2 &namdrg_bot ! BOTTOM friction (ln_drg_OFF =F) 3 3 !----------------------------------------------------------------------- 4 4 rn_Cd0 = 1.e-3 ! drag coefficient [-] -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_top
r10075 r13877 1 1 !----------------------------------------------------------------------- 2 &namdrg_top ! TOP friction (ln_ OFF =F & ln_isfcav=T)2 &namdrg_top ! TOP friction (ln_drg_OFF =F & ln_isfcav=T) 3 3 !----------------------------------------------------------------------- 4 4 rn_Cd0 = 1.e-3 ! drag coefficient [-] -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn
r11703 r13877 10 10 rn_ishlat = 2. ! lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 11 11 ln_landfast_L16 = .false. ! landfast: parameterization from Lemieux 2016 12 rn_ depfra= 0.125 ! fraction of ocean depth that ice must reach to initiate landfast12 rn_lf_depfra = 0.125 ! fraction of ocean depth that ice must reach to initiate landfast 13 13 ! recommended range: [0.1 ; 0.25] 14 rn_ icebfr = 15. ! maximum bottom stress per unit volume [N/m3]15 rn_lf relax= 1.e-5 ! relaxation time scale to reach static friction [s-1]16 rn_ tensile= 0.05 ! isotropic tensile strength [0-0.5??]14 rn_lf_bfr = 15. ! maximum bottom stress per unit volume [N/m3] 15 rn_lf_relax = 1.e-5 ! relaxation time scale to reach static friction [s-1] 16 rn_lf_tensile = 0.05 ! isotropic tensile strength [0-0.5??] 17 17 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn_rhg
r11025 r13877 9 9 rn_relast = 0.333 ! ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 10 10 ! advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 11 ln_rhg_chkcvg = .false. ! check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg) 11 12 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nameos
r10075 r13877 7 7 ! 8 8 ! ! S-EOS coefficients (ln_seos=T): 9 ! ! rd(T,S,Z)*r au0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS9 ! ! rd(T,S,Z)*rho0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS 10 10 rn_a0 = 1.6550e-1 ! thermal expension coefficient 11 11 rn_b0 = 7.6554e-1 ! saline expension coefficient -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namini
r11703 r13877 3 3 !------------------------------------------------------------------------------ 4 4 ln_iceini = .true. ! activate ice initialization (T) or not (F) 5 ln_iceini_file = .false. ! netcdf file provided for initialization (T) or not (F) 5 nn_iceini_file = 0 ! 0 = Initialise sea ice based on SSTs 6 ! 1 = Initialise sea ice from single category netcdf file 7 ! 2 = Initialise sea ice from multi category restart file 6 8 rn_thres_sst = 2.0 ! max temp. above Tfreeze with initial ice = (sst - tfreeze) 7 9 rn_hti_ini_n = 3.0 ! initial ice thickness (m), North … … 23 25 rn_hpd_ini_n = 0.05 ! initial pond depth (m), North 24 26 rn_hpd_ini_s = 0.05 ! " " South 25 ! -- for ln_iceini_file = T 27 rn_hld_ini_n = 0.0 ! initial pond lid depth (m), North 28 rn_hld_ini_s = 0.0 ! " " South 29 ! -- for nn_iceini_file = 1 26 30 sn_hti = 'Ice_initialization' , -12 ,'hti' , .false. , .true., 'yearly' , '' , '', '' 27 31 sn_hts = 'Ice_initialization' , -12 ,'hts' , .false. , .true., 'yearly' , '' , '', '' … … 34 38 sn_apd = 'NOT USED' , -12 ,'apd' , .false. , .true., 'yearly' , '' , '', '' 35 39 sn_hpd = 'NOT USED' , -12 ,'hpd' , .false. , .true., 'yearly' , '' , '', '' 40 sn_hld = 'NOT USED' , -12 ,'hld' , .false. , .true., 'yearly' , '' , '', '' 36 41 cn_dir='./' 37 42 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_blk
r12377 r13877 35 35 sn_tair = 't_10.15JUNE2009_fill' , 6. , 'T_10_MOD', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 36 36 sn_humi = 'q_10.15JUNE2009_fill' , 6. , 'Q_10_MOD', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 37 sn_hpgi = 'NONE' , 24. , 'uhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG' , ''38 sn_hpgj = 'NONE' , 24. , 'vhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG' , ''39 37 sn_prec = 'ncar_precip.15JUNE2009_fill', -1. , 'PRC_MOD1', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 40 38 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 41 39 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 40 sn_uoatm = 'NOT USED' , 6. , 'UOATM' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , 'Uoceatm', '' 41 sn_voatm = 'NOT USED' , 6. , 'VOATM' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , 'Voceatm', '' 42 sn_cc = 'NOT USED' , 24. , 'CC' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 43 sn_hpgi = 'NOT USED' , 24. , 'uhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG' , '' 44 sn_hpgj = 'NOT USED' , 24. , 'vhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG' , '' 42 45 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_cpl
r10075 r13877 2 2 &namsbc_cpl ! coupled ocean/atmosphere model ("key_oasis3") 3 3 !----------------------------------------------------------------------- 4 nn_cplmodel = 1 ! Maximum number of models to/from which NEMO is potentially sending/receiving data5 ln_usecplmask = .false. ! use a coupling mask file to merge data received from several models6 ! ! -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)7 nn_cats_cpl = 5 ! Number of sea ice categories over which coupling is to be carried out (if not 1)8 4 nn_cplmodel = 1 ! Maximum number of models to/from which NEMO is potentially sending/receiving data 5 ln_usecplmask = .false. ! use a coupling mask file to merge data received from several models 6 ! ! -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 7 ln_scale_ice_flux = .false. ! use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration) 8 nn_cats_cpl = 5 ! Number of sea ice categories over which coupling is to be carried out (if not 1) 9 9 !_____________!__________________________!____________!_____________!______________________!________! 10 10 ! ! description ! multiple ! vector ! vector ! vector ! -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd
r11025 r13877 6 6 ln_icedO = .true. ! activate ice growth in open-water (T) or not (F) 7 7 ln_icedS = .true. ! activate brine drainage (T) or not (F) 8 ! 9 ln_leadhfx = .true. ! heat in the leads is used to melt sea-ice before warming the ocean 8 10 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_pnd
r11536 r13877 2 2 &namthd_pnd ! Melt ponds 3 3 !------------------------------------------------------------------------------ 4 ln_pnd = .false. ! activate melt ponds or not 5 ln_pnd_H12 = .false. ! activate evolutive melt ponds (from Holland et al 2012) 6 ln_pnd_CST = .false. ! activate constant melt ponds 7 rn_apnd = 0.2 ! prescribed pond fraction, at Tsu=0 degC 8 rn_hpnd = 0.05 ! prescribed pond depth, at Tsu=0 degC 9 ln_pnd_alb = .false. ! melt ponds affect albedo or not 4 ln_pnd = .false. ! activate melt ponds or not 5 ln_pnd_LEV = .false. ! level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012) 6 rn_apnd_min = 0.15 ! minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ?? 7 rn_apnd_max = 0.85 ! maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ?? 8 ln_pnd_CST = .false. ! constant melt ponds 9 rn_apnd = 0.2 ! prescribed pond fraction, at Tsu=0 degC 10 rn_hpnd = 0.05 ! prescribed pond depth, at Tsu=0 degC 11 ln_pnd_lids = .true. ! frozen lids on top of the ponds (only for ln_pnd_LEV) 12 ln_pnd_alb = .true. ! effect of melt ponds on ice albedo 10 13 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_zdf
r11025 r13877 7 7 rn_cnd_s = 0.31 ! thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 8 8 ! Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 9 rn_kappa_i = 1.0 ! radiation attenuation coefficient in sea ice [1/m] 9 rn_kappa_i = 1.0 ! radiation attenuation coefficient in sea ice [1/m] 10 rn_kappa_s = 10.0 ! nn_qtrice = 0: radiation attenuation coefficient in snow [1/m] 11 rn_kappa_smlt = 7.0 ! nn_qtrice = 1: radiation attenuation coefficient in melting snow [1/m] 12 rn_kappa_sdry = 10.0 ! radiation attenuation coefficient in dry snow [1/m] 13 ln_zdf_chkcvg = .false. ! check convergence of heat diffusion scheme (output variable: tice_cvg) 10 14 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_gls
r9355 r13877 13 13 nn_z0_met = 2 ! Method for surface roughness computation (0/1/2/3) 14 14 ! ! =3 requires ln_wave=T 15 nn_z0_ice = 1 ! attenutaion of surface wave breaking under ice 16 ! ! = 0 no impact of ice cover 17 ! ! = 1 roughness uses rn_hsri and is weigthed by 1-TANH(10*fr_i) 18 ! ! = 2 roughness uses rn_hsri and is weighted by 1-fr_i 19 ! ! = 3 roughness uses rn_hsri and is weighted by 1-MIN(1,4*fr_i) 15 20 nn_bc_surf = 1 ! surface condition (0/1=Dir/Neum) 16 21 nn_bc_bot = 1 ! bottom condition (0/1=Dir/Neum) -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_tke
r10075 r13877 15 15 ln_mxl0 = .true. ! surface mixing length scale = F(wind stress) (T) or not (F) 16 16 rn_mxl0 = 0.04 ! surface buoyancy lenght scale minimum value 17 ln_drg = .false. ! top/bottom friction added as boundary condition of TKE18 17 ln_lc = .true. ! Langmuir cell parameterisation (Axell 2002) 19 18 rn_lc = 0.15 ! coef. associated to Langmuir cells … … 26 25 ! = 0 constant 10 m length scale 27 26 ! = 1 0.5m at the equator to 30m poleward of 40 degrees 28 rn_eice = 4 ! below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4 27 nn_eice = 1 ! attenutaion of langmuir & surface wave breaking under ice 28 ! ! = 0 no impact of ice cover on langmuir & surface wave breaking 29 ! ! = 1 weigthed by 1-TANH(10*fr_i) 30 ! ! = 2 weighted by 1-fr_i 31 ! ! = 3 weighted by 1-MIN(1,4*fr_i) 29 32 / -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice.F90
r13571 r13877 391 391 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vice !: ice volume variation [m/s] 392 392 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 393 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_aice !: ice conc. variation [s-1] 394 ! 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_mass !: advection of mass (kg/m2/s) 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_salt !: advection of salt (g/m2/s) 397 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_heat !: advection of heat (W/m2) 393 398 ! 394 399 !!---------------------------------------------------------------------- … … 493 498 ! * Ice diagnostics 494 499 ii = ii + 1 495 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 496 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 497 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), STAT=ierr(ii) ) 500 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 501 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 502 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), & 503 & diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 498 504 499 505 ! * Ice conservation -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice1d.F90
r13571 r13877 145 145 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sst_1d 146 146 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sss_1d 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: frq_m_1d 147 148 148 149 ! convergence check … … 225 226 ! 226 227 ii = ii + 1 227 ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) )228 ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , frq_m_1d(jpij) , STAT=ierr(ii) ) 228 229 ! 229 230 ii = ii + 1 -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icecor.F90
r13630 r13877 55 55 INTEGER :: ji, jj, jk, jl ! dummy loop indices 56 56 REAL(wp) :: zsal, zzc 57 REAL(wp), DIMENSION(jpi,jpj) :: zafx ! concentration trends diag58 57 !!---------------------------------------------------------------------- 59 58 ! controls … … 95 94 zsal = sv_i(ji,jj,jl) 96 95 sv_i(ji,jj,jl) = MIN( MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl) ) 97 sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 96 IF( kn /= 0 ) & ! no ice-ocean exchanges if kn=0 (for bdy for instance) otherwise conservation diags will fail 97 & sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc ! associated salt flux 98 98 END_2D 99 99 END DO 100 100 ENDIF 101 ! !-----------------------------------------------------102 CALL ice_var_zapsmall ! Zap small values !103 ! !-----------------------------------------------------104 101 102 IF( kn /= 0 ) THEN ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx) 103 ! otherwise conservation diags will fail 104 ! !----------------------------------------------------- 105 CALL ice_var_zapsmall ! Zap small values ! 106 ! !----------------------------------------------------- 107 ENDIF 105 108 ! !----------------------------------------------------- 106 109 IF( kn == 2 ) THEN ! Ice drift case: Corrections to avoid wrong values ! … … 119 122 #endif 120 123 ENDIF 121 122 ! !-----------------------------------------------------123 SELECT CASE( kn ) ! Diagnostics !124 ! !-----------------------------------------------------125 CASE( 1 ) !--- dyn trend diagnostics126 !127 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN128 diag_heat(:,:) = - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & ! W.m-2129 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice130 diag_sice(:,:) = SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi131 diag_vice(:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi132 diag_vsnw(:,:) = SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos133 ENDIF134 ! ! concentration tendency (dynamics)135 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN136 zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice137 CALL iom_put( 'afxdyn' , zafx )138 ENDIF139 !140 CASE( 2 ) !--- thermo trend diagnostics & ice aging141 !142 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! ice natural aging incrementation143 !144 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN145 diag_heat(:,:) = diag_heat(:,:) &146 & - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice &147 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice148 diag_sice(:,:) = diag_sice(:,:) &149 & + SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi150 diag_vice(:,:) = diag_vice(:,:) &151 & + SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi152 diag_vsnw(:,:) = diag_vsnw(:,:) &153 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos154 CALL iom_put ( 'hfxdhc' , diag_heat )155 ENDIF156 ! ! concentration tendency (total + thermo)157 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN158 zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice159 CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice )160 CALL iom_put( 'afxtot' , zafx )161 ENDIF162 !163 END SELECT164 124 ! 165 125 ! controls -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icectl.F90
r13571 r13877 43 43 PUBLIC ice_prt 44 44 PUBLIC ice_prt3D 45 PUBLIC ice_drift_wri 46 PUBLIC ice_drift_init 45 47 46 48 ! thresold rates for conservation … … 49 51 REAL(wp), PARAMETER :: zchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg) 50 52 REAL(wp), PARAMETER :: zchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg) 53 54 ! for drift outputs 55 CHARACTER(LEN=50) :: clname="icedrift_diagnostics.ascii" ! ascii filename 56 INTEGER :: numicedrift ! outfile unit 57 REAL(wp) :: rdiag_icemass, rdiag_icesalt, rdiag_iceheat 58 REAL(wp) :: rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat 51 59 52 60 !! * Substitutions … … 132 140 133 141 ! -- advection scheme is conservative? -- ! 134 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather)135 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) ! must be close to 0 (only for Prather)142 zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t ) 143 zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 136 144 137 145 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 156 164 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_amax 157 165 ! check if advection scheme is conservative 158 ! only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 159 ! so the formulation for conservation is different (and not coded) 160 ! it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 161 !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 162 ! & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 166 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 167 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 168 IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 169 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * rdt_ice 163 170 ENDIF 164 171 ! … … 186 193 ! water flux 187 194 ! -- mass diag -- ! 188 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 195 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 196 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) 189 197 190 198 ! -- salt diag -- ! 191 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t )199 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) 192 200 193 201 ! -- heat diag -- ! 194 ! clem: not the good formulation 195 !!zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr & 196 !! & ) * e1e2t ) 202 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 203 ! equivalent to this: 204 !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 205 !! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr & 206 !! & ) * e1e2t ) 197 207 198 208 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 204 214 IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 205 215 & WRITE(numout,*) cd_routine,' : violation salt cons. [g] = ',zdiag_salt * rDt_ice 206 !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 216 IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 217 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 207 218 ENDIF 208 219 ! … … 725 736 726 737 END SUBROUTINE ice_prt3D 738 739 740 SUBROUTINE ice_drift_wri( kt ) 741 !!------------------------------------------------------------------- 742 !! *** ROUTINE ice_drift_wri *** 743 !! 744 !! ** Purpose : conservation of mass, salt and heat 745 !! write the drift in a ascii file at each time step 746 !! and the total run drifts 747 !!------------------------------------------------------------------- 748 INTEGER, INTENT(in) :: kt ! ice time-step index 749 ! 750 INTEGER :: ji, jj 751 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 752 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass2D, zdiag_salt2D, zdiag_heat2D 753 !!------------------------------------------------------------------- 754 ! 755 IF( kt == nit000 .AND. lwp ) THEN 756 WRITE(numout,*) 757 WRITE(numout,*) 'ice_drift_wri: sea-ice drifts' 758 WRITE(numout,*) '~~~~~~~~~~~~~' 759 ENDIF 760 ! 761 ! 2D budgets (must be close to 0) 762 IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN 763 DO_2D( 1, 1, 1, 1 ) 764 zdiag_mass2D(ji,jj) = wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) & 765 & + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj) 766 zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj) 767 zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj) 768 END_2D 769 ! 770 ! write outputs 771 CALL iom_put( 'icedrift_mass', zdiag_mass2D ) 772 CALL iom_put( 'icedrift_salt', zdiag_salt2D ) 773 CALL iom_put( 'icedrift_heat', zdiag_heat2D ) 774 ENDIF 775 776 ! -- mass diag -- ! 777 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 778 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 779 zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 780 781 ! -- salt diag -- ! 782 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 783 zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 784 785 ! -- heat diag -- ! 786 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 787 zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t ) 788 789 ! ! write out to file 790 IF( lwp ) THEN 791 ! check global drift (must be close to 0) 792 WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift [kg]', zdiag_mass 793 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift [kg]', zdiag_salt 794 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift [W] ', zdiag_heat 795 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 796 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'mass drift adv [kg]', zdiag_adv_mass 797 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift adv [kg]', zdiag_adv_salt 798 WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift adv [W] ', zdiag_adv_heat 799 ENDIF 800 ! ! drifts 801 rdiag_icemass = rdiag_icemass + zdiag_mass 802 rdiag_icesalt = rdiag_icesalt + zdiag_salt 803 rdiag_iceheat = rdiag_iceheat + zdiag_heat 804 rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass 805 rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt 806 rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat 807 ! 808 ! ! output drifts and close ascii file 809 IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN 810 ! to ascii file 811 WRITE(numicedrift,*) '******************************************' 812 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift [kg]', rdiag_icemass 813 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass 814 WRITE(numicedrift,*) '******************************************' 815 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift [kg]', rdiag_icesalt 816 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt 817 WRITE(numicedrift,*) '******************************************' 818 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift [W] ', rdiag_iceheat 819 WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat 820 CLOSE( numicedrift ) 821 ! 822 ! to ocean output 823 WRITE(numout,*) 824 WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run ' 825 WRITE(numout,*) '~~~~~~~~~~~~~' 826 ! check global drift (must be close to 0) 827 WRITE(numout,*) ' sea-ice mass drift [kg] = ', rdiag_icemass 828 WRITE(numout,*) ' sea-ice salt drift [kg] = ', rdiag_icesalt 829 WRITE(numout,*) ' sea-ice heat drift [W] = ', rdiag_iceheat 830 ! check drift from advection scheme (can be /=0 with bdy but not sure why) 831 WRITE(numout,*) ' sea-ice mass drift adv [kg] = ', rdiag_adv_icemass 832 WRITE(numout,*) ' sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt 833 WRITE(numout,*) ' sea-ice heat drift adv [W] = ', rdiag_adv_iceheat 834 ENDIF 835 ! 836 END SUBROUTINE ice_drift_wri 837 838 SUBROUTINE ice_drift_init 839 !!---------------------------------------------------------------------- 840 !! *** ROUTINE ice_drift_init *** 841 !! 842 !! ** Purpose : create output file, initialise arrays 843 !!---------------------------------------------------------------------- 844 ! 845 IF( .NOT.ln_icediachk ) RETURN ! exit 846 ! 847 IF(lwp) THEN 848 WRITE(numout,*) 849 WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file' 850 WRITE(numout,*) '~~~~~~~~~~~~~' 851 WRITE(numout,*) 852 ! 853 ! create output ascii file 854 CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea ) 855 WRITE(numicedrift,*) 'Timestep Drifts' 856 WRITE(numicedrift,*) '******************************************' 857 ENDIF 858 ! 859 rdiag_icemass = 0._wp 860 rdiag_icesalt = 0._wp 861 rdiag_iceheat = 0._wp 862 rdiag_adv_icemass = 0._wp 863 rdiag_adv_icesalt = 0._wp 864 rdiag_adv_iceheat = 0._wp 865 ! 866 END SUBROUTINE ice_drift_init 727 867 728 868 #else -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_pra.F90
r13807 r13877 88 88 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 89 89 INTEGER :: icycle ! number of sub-timestep for the advection 90 REAL(wp) :: zdt 90 REAL(wp) :: zdt, z1_dt ! - - 91 91 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 92 92 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 … … 100 100 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 101 101 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei 102 !! diagnostics 103 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 102 104 !!---------------------------------------------------------------------- 103 105 ! … … 109 111 ELSEWHERE ; zs_i(:,:,:) = 0._wp 110 112 END WHERE 111 DO jl = 1, jpl 112 DO_2D( 0, 0, 0, 0 ) 113 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 114 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 115 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 116 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 117 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 118 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 119 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 120 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 121 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 122 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 123 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 124 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 125 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 126 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 127 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 128 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 129 END_2D 130 END DO 113 CALL icemax3D( ph_i , zhi_max ) 114 CALL icemax3D( ph_s , zhs_max ) 115 CALL icemax3D( ph_ip, zhip_max) 116 CALL icemax3D( zs_i , zsi_max ) 131 117 #if defined key_mpi3 132 118 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) … … 145 131 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 146 132 END WHERE 147 END DO 148 DO jl = 1, jpl 149 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 150 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 151 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 152 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 153 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 154 END_3D 155 END DO 156 DO jl = 1, jpl 157 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 158 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 159 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 160 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 161 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 162 END_3D 163 END DO 133 END DO 134 CALL icemax4D( ze_i , zei_max ) 135 CALL icemax4D( ze_s , zes_max ) 164 136 #if defined key_mpi3 165 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1. )166 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1. )137 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 138 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 167 139 #else 168 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. )169 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. )140 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 141 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 170 142 #endif 171 143 ! … … 184 156 ENDIF 185 157 zdt = rDt_ice / REAL(icycle) 158 z1_dt = 1._wp / zdt 186 159 187 160 ! --- transport --- ! … … 190 163 191 164 DO jt = 1, icycle 165 166 ! diagnostics 167 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 168 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 169 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 170 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 192 171 193 172 ! record at_i before advection (for open water) … … 288 267 CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 289 268 ENDIF 290 ENDIF 291 ! 269 ENDIF 270 ! 271 ENDIF 272 273 ! --- Lateral boundary conditions --- ! 274 ! caution: for gradients (sx and sy) the sign changes 275 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp & ! ice volume 276 & , sxxice, 'T', 1._wp, syyice, 'T', 1._wp, sxyice, 'T', 1._wp & 277 & , z0snw , 'T', 1._wp, sxsn , 'T', -1._wp, sysn , 'T', -1._wp & ! snw volume 278 & , sxxsn , 'T', 1._wp, syysn , 'T', 1._wp, sxysn , 'T', 1._wp ) 279 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp & ! ice salinity 280 & , sxxsal, 'T', 1._wp, syysal, 'T', 1._wp, sxysal, 'T', 1._wp & 281 & , z0ai , 'T', 1._wp, sxa , 'T', -1._wp, sya , 'T', -1._wp & ! ice concentration 282 & , sxxa , 'T', 1._wp, syya , 'T', 1._wp, sxya , 'T', 1._wp ) 283 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp & ! ice age 284 & , sxxage, 'T', 1._wp, syyage, 'T', 1._wp, sxyage, 'T', 1._wp ) 285 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es , 'T', 1._wp, sxc0 , 'T', -1._wp, syc0 , 'T', -1._wp & ! snw enthalpy 286 & , sxxc0 , 'T', 1._wp, syyc0 , 'T', 1._wp, sxyc0 , 'T', 1._wp ) 287 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy 288 & , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp ) 289 IF ( ln_pnd_LEV ) THEN 290 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction 291 & , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp & 292 & , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp & ! melt pond volume 293 & , sxxvp, 'T', 1._wp, syyvp, 'T', 1._wp, sxyvp, 'T', 1._wp ) 294 IF ( ln_pnd_lids ) THEN 295 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp & ! melt pond lid volume 296 & , sxxvl,'T', 1._wp, syyvl,'T', 1._wp, sxyvl,'T', 1._wp ) 297 ENDIF 292 298 ENDIF 293 299 … … 325 331 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp ) 326 332 #endif 333 ! 334 ! --- diagnostics --- ! 335 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 336 & - zdiag_adv_mass(:,:) ) * z1_dt 337 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 338 & - zdiag_adv_salt(:,:) ) * z1_dt 339 diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 340 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 341 & - zdiag_adv_heat(:,:) ) * z1_dt 327 342 ! 328 343 ! --- Ensure non-negative fields --- ! … … 363 378 !! 364 379 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 380 INTEGER :: jj0 ! dummy loop indices 365 381 REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars 366 382 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 367 383 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 384 REAL(wp) :: zpsm, zps0 385 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy 368 386 REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace 369 387 REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - - 370 388 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 371 389 !----------------------------------------------------------------------- 390 ! in order to avoid lbc_lnk (communications): 391 ! jj loop must be 1:jpj if adv_x is called first 392 ! and 2:jpj-1 if adv_x is called second 393 jj0 = NINT(pcrh) 372 394 ! 373 395 jcat = SIZE( ps0 , 3 ) ! size of input arrays … … 376 398 ! 377 399 ! Limitation of moments. 378 DO_2D( 0, 0, 1, 1 ) 379 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 380 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 381 ! 382 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 383 zs1max = 1.5 * zslpmax 384 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 385 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 386 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 387 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 388 389 ps0 (ji,jj,jl) = zslpmax 390 psx (ji,jj,jl) = zs1new * rswitch 391 psxx(ji,jj,jl) = zs2new * rswitch 392 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 393 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 394 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 395 END_2D 396 397 ! Calculate fluxes and moments between boxes i<-->i+1 398 DO_2D( 0, 0, 1, 1 ) ! Flux from i to i+1 WHEN u GT 0 399 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 400 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 401 zalfq = zalf * zalf 402 zalf1 = 1.0 - zalf 403 zalf1q = zalf1 * zalf1 404 ! 405 zfm (ji,jj) = zalf * psm (ji,jj,jl) 406 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 407 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 408 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 409 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 410 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 411 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 412 413 ! Readjust moments remaining in the box. 414 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 415 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 416 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 417 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 418 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 419 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 420 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 421 END_2D 422 423 DO_2D( 0, 0, 1, 0 ) ! Flux from i+1 to i when u LT 0. 424 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 425 zalg (ji,jj) = zalf 426 zalfq = zalf * zalf 427 zalf1 = 1.0 - zalf 428 zalg1 (ji,jj) = zalf1 429 zalf1q = zalf1 * zalf1 430 zalg1q(ji,jj) = zalf1q 431 ! 432 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 433 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 434 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 435 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 436 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 437 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 438 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 439 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 440 END_2D 441 442 DO_2D( 0, 0, 0, 0 ) ! Readjust moments remaining in the box. 443 zbt = zbet(ji-1,jj) 444 zbt1 = 1.0 - zbet(ji-1,jj) 445 ! 446 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 447 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 448 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 449 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 450 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 451 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 452 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 453 END_2D 454 455 ! Put the temporary moments into appropriate neighboring boxes. 456 DO_2D( 0, 0, 0, 0 ) ! Flux from i to i+1 IF u GT 0. 457 zbt = zbet(ji-1,jj) 458 zbt1 = 1.0 - zbet(ji-1,jj) 459 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 460 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 461 zalf1 = 1.0 - zalf 462 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 463 ! 464 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 465 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 466 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 467 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 468 & + zbt1 * psxx(ji,jj,jl) 469 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 470 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 471 & + zbt1 * psxy(ji,jj,jl) 472 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 473 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 474 END_2D 475 476 DO_2D( 0, 0, 0, 0 ) ! Flux from i+1 to i IF u LT 0. 477 zbt = zbet(ji,jj) 478 zbt1 = 1.0 - zbet(ji,jj) 479 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 480 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 481 zalf1 = 1.0 - zalf 482 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 483 ! 484 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 485 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 486 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 487 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 488 & + ( zalf1 - zalf ) * ztemp ) ) 489 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 490 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 491 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 492 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 493 END_2D 494 400 DO jj = Njs0 - jj0, Nje0 + jj0 401 402 DO ji = Nis0 - 1, Nie0 + 1 403 404 zpsm = psm (ji,jj,jl) ! optimization 405 zps0 = ps0 (ji,jj,jl) 406 zpsx = psx (ji,jj,jl) 407 zpsxx = psxx(ji,jj,jl) 408 zpsy = psy (ji,jj,jl) 409 zpsyy = psyy(ji,jj,jl) 410 zpsxy = psxy(ji,jj,jl) 411 412 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 413 zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 414 ! 415 zslpmax = MAX( 0._wp, zps0 ) 416 zs1max = 1.5 * zslpmax 417 zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) ) 418 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 419 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 420 421 zps0 = zslpmax 422 zpsx = zs1new * rswitch 423 zpsxx = zs2new * rswitch 424 zpsy = zpsy * rswitch 425 zpsyy = zpsyy * rswitch 426 zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 427 428 ! Calculate fluxes and moments between boxes i<-->i+1 429 ! ! Flux from i to i+1 WHEN u GT 0 430 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 431 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 432 zalfq = zalf * zalf 433 zalf1 = 1.0 - zalf 434 zalf1q = zalf1 * zalf1 435 ! 436 zfm (ji,jj) = zalf * zpsm 437 zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 438 zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx ) 439 zfxx(ji,jj) = zalf * zpsxx * zalfq 440 zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy ) 441 zfxy(ji,jj) = zalfq * zpsxy 442 zfyy(ji,jj) = zalf * zpsyy 443 444 ! ! Readjust moments remaining in the box. 445 zpsm = zpsm - zfm(ji,jj) 446 zps0 = zps0 - zf0(ji,jj) 447 zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 448 zpsxx = zalf1 * zalf1q * zpsxx 449 zpsy = zpsy - zfy (ji,jj) 450 zpsyy = zpsyy - zfyy(ji,jj) 451 zpsxy = zalf1q * zpsxy 452 ! 453 psm (ji,jj,jl) = zpsm ! optimization 454 ps0 (ji,jj,jl) = zps0 455 psx (ji,jj,jl) = zpsx 456 psxx(ji,jj,jl) = zpsxx 457 psy (ji,jj,jl) = zpsy 458 psyy(ji,jj,jl) = zpsyy 459 psxy(ji,jj,jl) = zpsxy 460 ! 461 END DO 462 463 DO ji = Nis0 - 1, Nie0 464 ! ! Flux from i+1 to i when u LT 0. 465 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 466 zalg (ji,jj) = zalf 467 zalfq = zalf * zalf 468 zalf1 = 1.0 - zalf 469 zalg1 (ji,jj) = zalf1 470 zalf1q = zalf1 * zalf1 471 zalg1q(ji,jj) = zalf1q 472 ! 473 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 474 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 475 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 476 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 477 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 478 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 479 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 480 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 481 END DO 482 483 DO ji = Nis0, Nie0 484 ! 485 zpsm = psm (ji,jj,jl) ! optimization 486 zps0 = ps0 (ji,jj,jl) 487 zpsx = psx (ji,jj,jl) 488 zpsxx = psxx(ji,jj,jl) 489 zpsy = psy (ji,jj,jl) 490 zpsyy = psyy(ji,jj,jl) 491 zpsxy = psxy(ji,jj,jl) 492 ! ! Readjust moments remaining in the box. 493 zbt = zbet(ji-1,jj) 494 zbt1 = 1.0 - zbet(ji-1,jj) 495 ! 496 zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 497 zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 498 zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 499 zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 500 zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) ) 501 zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 502 zpsxy = zalg1q(ji-1,jj) * zpsxy 503 504 ! Put the temporary moments into appropriate neighboring boxes. 505 ! ! Flux from i to i+1 IF u GT 0. 506 zbt = zbet(ji-1,jj) 507 zbt1 = 1.0 - zbet(ji-1,jj) 508 zpsm = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 509 zalf = zbt * zfm(ji-1,jj) / zpsm 510 zalf1 = 1.0 - zalf 511 ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 512 ! 513 zps0 = zbt * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 514 zpsx = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 515 zpsxx = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx & 516 & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 517 & + zbt1 * zpsxx 518 zpsxy = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy & 519 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * zpsy ) ) & 520 & + zbt1 * zpsxy 521 zpsy = zbt * ( zpsy + zfy (ji-1,jj) ) + zbt1 * zpsy 522 zpsyy = zbt * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 523 524 ! ! Flux from i+1 to i IF u LT 0. 525 zbt = zbet(ji,jj) 526 zbt1 = 1.0 - zbet(ji,jj) 527 zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 528 zalf = zbt1 * zfm(ji,jj) / zpsm 529 zalf1 = 1.0 - zalf 530 ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 531 ! 532 zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) 533 zpsx = zbt * zpsx + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 534 zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 535 & + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) ) & 536 & + ( zalf1 - zalf ) * ztemp ) ) 537 zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy & 538 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 539 zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) ) 540 zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 541 ! 542 psm (ji,jj,jl) = zpsm ! optimization 543 ps0 (ji,jj,jl) = zps0 544 psx (ji,jj,jl) = zpsx 545 psxx(ji,jj,jl) = zpsxx 546 psy (ji,jj,jl) = zpsy 547 psyy(ji,jj,jl) = zpsyy 548 psxy(ji,jj,jl) = zpsxy 549 END DO 550 ! 551 END DO 552 ! 495 553 END DO 496 497 !-- Lateral boundary conditions 498 #if defined key_mpi3 499 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp & 500 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes 501 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp ) 502 #else 503 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp & 504 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes 505 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp ) 506 #endif 507 ! 554 ! 508 555 END SUBROUTINE adv_x 509 556 … … 526 573 !! 527 574 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 575 INTEGER :: ji0 ! dummy loop indices 528 576 REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars 529 577 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 530 578 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 579 REAL(wp) :: zpsm, zps0 580 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy 531 581 REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace 532 582 REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - 533 583 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 534 584 !--------------------------------------------------------------------- 585 ! in order to avoid lbc_lnk (communications): 586 ! ji loop must be 1:jpi if adv_y is called first 587 ! and 2:jpi-1 if adv_y is called second 588 ji0 = NINT(pcrh) 535 589 ! 536 590 jcat = SIZE( ps0 , 3 ) ! size of input arrays … … 539 593 ! 540 594 ! Limitation of moments. 541 DO_2D( 1, 1, 0, 0 ) 542 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 543 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 544 ! 545 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 595 DO_2D( 1, 1, ji0, ji0 ) 596 ! 597 zpsm = psm (ji,jj,jl) ! optimization 598 zps0 = ps0 (ji,jj,jl) 599 zpsx = psx (ji,jj,jl) 600 zpsxx = psxx(ji,jj,jl) 601 zpsy = psy (ji,jj,jl) 602 zpsyy = psyy(ji,jj,jl) 603 zpsxy = psxy(ji,jj,jl) 604 ! 605 ! Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 606 zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 607 ! 608 zslpmax = MAX( 0._wp, zps0 ) 546 609 zs1max = 1.5 * zslpmax 547 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 548 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 549 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 610 zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) ) 611 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 550 612 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 551 613 ! 552 ps0 (ji,jj,jl) = zslpmax 553 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 554 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 555 psy (ji,jj,jl) = zs1new * rswitch 556 psyy(ji,jj,jl) = zs2new * rswitch 557 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 558 END_2D 559 560 ! Calculate fluxes and moments between boxes j<-->j+1 561 DO_2D( 1, 1, 0, 0 ) ! Flux from j to j+1 WHEN v GT 0 614 zps0 = zslpmax 615 zpsx = zpsx * rswitch 616 zpsxx = zpsxx * rswitch 617 zpsy = zs1new * rswitch 618 zpsyy = zs2new * rswitch 619 zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 620 621 ! Calculate fluxes and moments between boxes j<-->j+1 622 ! ! Flux from j to j+1 WHEN v GT 0 562 623 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 563 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)624 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 564 625 zalfq = zalf * zalf 565 626 zalf1 = 1.0 - zalf 566 627 zalf1q = zalf1 * zalf1 567 628 ! 568 zfm (ji,jj) = zalf * psm(ji,jj,jl) 569 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 570 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 571 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 572 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 573 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 574 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 575 ! 576 ! Readjust moments remaining in the box. 577 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 578 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 579 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 580 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 581 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 582 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 583 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 629 zfm (ji,jj) = zalf * zpsm 630 zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1-zalf) * zpsyy ) ) 631 zfy (ji,jj) = zalfq *( zpsy + 3.0*zalf1*zpsyy ) 632 zfyy(ji,jj) = zalf * zalfq * zpsyy 633 zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy ) 634 zfxy(ji,jj) = zalfq * zpsxy 635 zfxx(ji,jj) = zalf * zpsxx 636 ! 637 ! ! Readjust moments remaining in the box. 638 zpsm = zpsm - zfm(ji,jj) 639 zps0 = zps0 - zf0(ji,jj) 640 zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 641 zpsyy = zalf1 * zalf1q * zpsyy 642 zpsx = zpsx - zfx(ji,jj) 643 zpsxx = zpsxx - zfxx(ji,jj) 644 zpsxy = zalf1q * zpsxy 645 ! 646 psm (ji,jj,jl) = zpsm ! optimization 647 ps0 (ji,jj,jl) = zps0 648 psx (ji,jj,jl) = zpsx 649 psxx(ji,jj,jl) = zpsxx 650 psy (ji,jj,jl) = zpsy 651 psyy(ji,jj,jl) = zpsyy 652 psxy(ji,jj,jl) = zpsxy 584 653 END_2D 585 654 ! 586 DO_2D( 1, 0, 0, 0 ) ! Flux from j+1 to j when v LT 0. 655 DO_2D( 1, 0, ji0, ji0 ) 656 ! ! Flux from j+1 to j when v LT 0. 587 657 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 588 658 zalg (ji,jj) = zalf … … 603 673 END_2D 604 674 605 ! Readjust moments remaining in the box.606 DO_2D( 0, 0, 0, 0 )675 DO_2D( 0, 0, ji0, ji0 ) 676 ! ! Readjust moments remaining in the box. 607 677 zbt = zbet(ji,jj-1) 608 678 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 609 679 ! 610 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 611 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 612 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 613 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 614 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 615 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 616 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 680 zpsm = psm (ji,jj,jl) ! optimization 681 zps0 = ps0 (ji,jj,jl) 682 zpsx = psx (ji,jj,jl) 683 zpsxx = psxx(ji,jj,jl) 684 zpsy = psy (ji,jj,jl) 685 zpsyy = psyy(ji,jj,jl) 686 zpsxy = psxy(ji,jj,jl) 687 ! 688 zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 689 zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 690 zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 691 zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 692 zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) ) 693 zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 694 zpsxy = zalg1q(ji,jj-1) * zpsxy 695 696 ! Put the temporary moments into appropriate neighboring boxes. 697 ! ! Flux from j to j+1 IF v GT 0. 698 zbt = zbet(ji,jj-1) 699 zbt1 = 1.0 - zbet(ji,jj-1) 700 zpsm = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm 701 zalf = zbt * zfm(ji,jj-1) / zpsm 702 zalf1 = 1.0 - zalf 703 ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 704 ! 705 zps0 = zbt * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 706 zpsy = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp ) & 707 & + zbt1 * zpsy 708 zpsyy = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy & 709 & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 710 & + zbt1 * zpsyy 711 zpsxy = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy & 712 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) & 713 & + zbt1 * zpsxy 714 zpsx = zbt * ( zpsx + zfx (ji,jj-1) ) + zbt1 * zpsx 715 zpsxx = zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 716 717 ! ! Flux from j+1 to j IF v LT 0. 718 zbt = zbet(ji,jj) 719 zbt1 = 1.0 - zbet(ji,jj) 720 zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 721 zalf = zbt1 * zfm(ji,jj) / zpsm 722 zalf1 = 1.0 - zalf 723 ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 724 ! 725 zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) 726 zpsy = zbt * zpsy + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 727 zpsyy = zbt * zpsyy + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 728 & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) & 729 & + ( zalf1 - zalf ) * ztemp ) ) 730 zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy & 731 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 732 zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) ) 733 zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 734 ! 735 psm (ji,jj,jl) = zpsm ! optimization 736 ps0 (ji,jj,jl) = zps0 737 psx (ji,jj,jl) = zpsx 738 psxx(ji,jj,jl) = zpsxx 739 psy (ji,jj,jl) = zpsy 740 psyy(ji,jj,jl) = zpsyy 741 psxy(ji,jj,jl) = zpsxy 617 742 END_2D 618 619 ! Put the temporary moments into appropriate neighboring boxes. 620 DO_2D( 0, 0, 0, 0 ) ! Flux from j to j+1 IF v GT 0. 621 zbt = zbet(ji,jj-1) 622 zbt1 = 1.0 - zbet(ji,jj-1) 623 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 624 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 625 zalf1 = 1.0 - zalf 626 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 627 ! 628 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 629 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 630 & + zbt1 * psy(ji,jj,jl) 631 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 632 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 633 & + zbt1 * psyy(ji,jj,jl) 634 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 635 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 636 & + zbt1 * psxy(ji,jj,jl) 637 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 638 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 639 END_2D 640 641 DO_2D( 0, 0, 0, 0 ) ! Flux from j+1 to j IF v LT 0. 642 zbt = zbet(ji,jj) 643 zbt1 = 1.0 - zbet(ji,jj) 644 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 645 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 646 zalf1 = 1.0 - zalf 647 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 648 ! 649 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 650 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 651 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 652 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 653 & + ( zalf1 - zalf ) * ztemp ) ) 654 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 655 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 656 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 657 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 658 END_2D 659 743 ! 660 744 END DO 661 662 !-- Lateral boundary conditions663 #if defined key_mpi3664 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp &665 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes666 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp )667 #else668 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp &669 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes670 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp )671 #endif672 745 ! 673 746 END SUBROUTINE adv_y … … 897 970 ! 898 971 ! ! ice thickness 899 CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice )900 CALL iom_get( numrir, jpdom_auto, 'syice' , syice )972 CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp ) 973 CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp ) 901 974 CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 902 975 CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 903 976 CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 904 977 ! ! snow thickness 905 CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn )906 CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn )978 CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn , psgn = -1._wp ) 979 CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn , psgn = -1._wp ) 907 980 CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn ) 908 981 CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn ) 909 982 CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn ) 910 983 ! ! ice concentration 911 CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa )912 CALL iom_get( numrir, jpdom_auto, 'sya' , sya )984 CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa , psgn = -1._wp ) 985 CALL iom_get( numrir, jpdom_auto, 'sya' , sya , psgn = -1._wp ) 913 986 CALL iom_get( numrir, jpdom_auto, 'sxxa' , sxxa ) 914 987 CALL iom_get( numrir, jpdom_auto, 'syya' , syya ) 915 988 CALL iom_get( numrir, jpdom_auto, 'sxya' , sxya ) 916 989 ! ! ice salinity 917 CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal )918 CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal )990 CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp ) 991 CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp ) 919 992 CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 920 993 CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 921 994 CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 922 995 ! ! ice age 923 CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage )924 CALL iom_get( numrir, jpdom_auto, 'syage' , syage )996 CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp ) 997 CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp ) 925 998 CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 926 999 CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) … … 929 1002 DO jk = 1, nlay_s 930 1003 WRITE(zchar1,'(I2.2)') jk 931 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)932 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:)1004 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxc0 (:,:,jk,:) = z3d(:,:,:) 1005 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; syc0 (:,:,jk,:) = z3d(:,:,:) 933 1006 znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:) 934 1007 znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:) … … 938 1011 DO jk = 1, nlay_i 939 1012 WRITE(zchar1,'(I2.2)') jk 940 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:)941 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:)1013 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxe (:,:,jk,:) = z3d(:,:,:) 1014 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sye (:,:,jk,:) = z3d(:,:,:) 942 1015 znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:) 943 1016 znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:) … … 946 1019 ! 947 1020 IF( ln_pnd_LEV ) THEN ! melt pond fraction 948 IF( iom_varid( numr or, 'sxap', ldstop = .FALSE. ) > 0 ) THEN949 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap )950 CALL iom_get( numrir, jpdom_auto, 'syap' , syap )1021 IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 1022 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) 1023 CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp ) 951 1024 CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 952 1025 CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 953 1026 CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 954 1027 ! ! melt pond volume 955 CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp )956 CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp )1028 CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp ) 1029 CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp ) 957 1030 CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 958 1031 CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) … … 964 1037 ! 965 1038 IF ( ln_pnd_lids ) THEN ! melt pond lid volume 966 IF( iom_varid( numr or, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN967 CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl )968 CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl )1039 IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 1040 CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp ) 1041 CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp ) 969 1042 CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 970 1043 CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) … … 1081 1154 END SUBROUTINE adv_pra_rst 1082 1155 1156 SUBROUTINE icemax3D( pice , pmax ) 1157 !!--------------------------------------------------------------------- 1158 !! *** ROUTINE icemax3D *** 1159 !! ** Purpose : compute the max of the 9 points around 1160 !!---------------------------------------------------------------------- 1161 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1162 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1163 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1164 INTEGER :: ji, jj, jl ! dummy loop indices 1165 !!---------------------------------------------------------------------- 1166 DO jl = 1, jpl 1167 DO jj = Njs0-1, Nje0+1 1168 DO ji = Nis0, Nie0 1169 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1170 END DO 1171 END DO 1172 DO jj = Njs0, Nje0 1173 DO ji = Nis0, Nie0 1174 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1175 END DO 1176 END DO 1177 END DO 1178 END SUBROUTINE icemax3D 1179 1180 SUBROUTINE icemax4D( pice , pmax ) 1181 !!--------------------------------------------------------------------- 1182 !! *** ROUTINE icemax4D *** 1183 !! ** Purpose : compute the max of the 9 points around 1184 !!---------------------------------------------------------------------- 1185 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1186 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1187 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1188 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1189 !!---------------------------------------------------------------------- 1190 jlay = SIZE( pice , 3 ) ! size of input arrays 1191 DO jl = 1, jpl 1192 DO jk = 1, jlay 1193 DO jj = Njs0-1, Nje0+1 1194 DO ji = Nis0, Nie0 1195 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1196 END DO 1197 END DO 1198 DO jj = Njs0, Nje0 1199 DO ji = Nis0, Nie0 1200 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1201 END DO 1202 END DO 1203 END DO 1204 END DO 1205 END SUBROUTINE icemax4D 1206 1083 1207 #else 1084 1208 !!---------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_umx.F90
r13807 r13877 92 92 INTEGER :: icycle ! number of sub-timestep for the advection 93 93 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 94 REAL(wp) :: zdt, z vi_cen94 REAL(wp) :: zdt, z1_dt, zvi_cen 95 95 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 96 96 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box … … 104 104 ! 105 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 106 !! diagnostics 107 REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat 106 108 !!---------------------------------------------------------------------- 107 109 ! … … 113 115 ELSEWHERE ; zs_i(:,:,:) = 0._wp 114 116 END WHERE 115 DO jl = 1, jpl 116 DO_2D( 0, 0, 0, 0 ) 117 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 118 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 119 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 120 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 121 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 122 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 123 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 124 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 125 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 126 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 127 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 128 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 129 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 130 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 131 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 132 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 133 END_2D 134 END DO 117 CALL icemax3D( ph_i , zhi_max ) 118 CALL icemax3D( ph_s , zhs_max ) 119 CALL icemax3D( ph_ip, zhip_max) 120 CALL icemax3D( zs_i , zsi_max ) 135 121 #if defined key_mpi3 136 122 CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) … … 149 135 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 150 136 END WHERE 151 END DO 152 DO jl = 1, jpl 153 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 154 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 155 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 156 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 157 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 158 END_3D 159 END DO 160 DO jl = 1, jpl 161 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 162 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 163 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 164 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 165 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 166 END_3D 167 END DO 168 #if defined key_mpi3 169 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1. ) 170 CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1. ) 171 #else 172 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 173 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 174 #endif 137 END DO 138 CALL icemax4D( ze_i , zei_max ) 139 CALL icemax4D( ze_s , zes_max ) 140 CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp ) 141 CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp ) 175 142 ! 176 143 ! … … 188 155 ENDIF 189 156 zdt = rDt_ice / REAL(icycle) 157 z1_dt = 1._wp / zdt 190 158 191 159 ! --- transport --- ! … … 216 184 !---------------! 217 185 DO jt = 1, icycle 186 187 ! diagnostics 188 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 189 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 190 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 191 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) 218 192 219 193 ! record at_i before advection (for open water) … … 386 360 ENDIF 387 361 ENDIF 362 363 ! --- Lateral boundary conditions --- ! 364 IF ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 365 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 366 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 367 ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 368 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 369 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 370 ELSE 371 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp ) 372 ENDIF 373 CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 374 CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 388 375 ! 389 376 !== Open water area ==! … … 399 386 #endif 400 387 ! 388 ! --- diagnostics --- ! 389 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 390 & - zdiag_adv_mass(:,:) ) * z1_dt 391 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 392 & - zdiag_adv_salt(:,:) ) * z1_dt 393 diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 394 & - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) & 395 & - zdiag_adv_heat(:,:) ) * z1_dt 401 396 ! 402 397 ! --- Ensure non-negative fields and in-bound thicknesses --- ! … … 458 453 !! work on H (and not V). It is partly related to the multi-category approach 459 454 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 460 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 461 !! since sv_i and e_i are still good. 455 !! concentration is small). We also limit S and T. 462 456 !!---------------------------------------------------------------------- 463 457 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) … … 503 497 IF( pamsk == 0._wp ) THEN 504 498 DO jl = 1, jpl 505 DO_2D( 1, 0, 1, 0 )499 DO_2D( 0, 0, 1, 0 ) 506 500 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 507 501 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) … … 512 506 ENDIF 513 507 ! 508 END_2D 509 DO_2D( 1, 0, 0, 0 ) 514 510 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 515 511 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) … … 550 546 IF( PRESENT( pua_ho ) ) THEN 551 547 DO jl = 1, jpl 552 DO_2D( 1, 0, 1, 0 ) 553 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 554 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 548 DO_2D( 0, 0, 1, 0 ) 549 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 550 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 551 END_2D 552 DO_2D( 1, 0, 0, 0 ) 553 pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 554 pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 555 555 END_2D 556 556 END DO … … 566 566 END_2D 567 567 END DO 568 #if defined key_mpi3569 CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', ptc, 'T', 1.0_wp )570 #else571 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1.0_wp )572 #endif573 568 ! 574 569 END SUBROUTINE adv_umx … … 609 604 ! 610 605 DO jl = 1, jpl !-- flux in x-direction 611 DO_2D( 1, 0, 1, 0 )606 DO_2D( 1, 1, 1, 0 ) 612 607 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 613 608 END_2D … … 615 610 ! 616 611 DO jl = 1, jpl !-- first guess of tracer from u-flux 617 DO_2D( 0, 0, 0, 0 )612 DO_2D( 1, 1, 0, 0 ) 618 613 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 619 614 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 622 617 END_2D 623 618 END DO 624 #if defined key_mpi3625 CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )626 #else627 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )628 #endif629 619 ! 630 620 DO jl = 1, jpl !-- flux in y-direction 631 DO_2D( 1, 0, 1, 0 )621 DO_2D( 1, 0, 0, 0 ) 632 622 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 633 623 END_2D … … 637 627 ! 638 628 DO jl = 1, jpl !-- flux in y-direction 639 DO_2D( 1, 0, 1, 0)629 DO_2D( 1, 0, 1, 1 ) 640 630 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 641 631 END_2D … … 643 633 ! 644 634 DO jl = 1, jpl !-- first guess of tracer from v-flux 645 DO_2D( 0, 0, 0, 0)635 DO_2D( 0, 0, 1, 1 ) 646 636 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 647 637 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 650 640 END_2D 651 641 END DO 652 #if defined key_mpi3653 CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )654 #else655 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )656 #endif657 642 ! 658 643 DO jl = 1, jpl !-- flux in x-direction 659 DO_2D( 1, 0, 1, 0 )644 DO_2D( 0, 0, 1, 0 ) 660 645 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 661 646 END_2D … … 710 695 ! 711 696 DO jl = 1, jpl 712 DO_2D( 1, 0, 1, 0 )697 DO_2D( 1, 1, 1, 0 ) 713 698 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 699 END_2D 700 DO_2D( 1, 0, 1, 1 ) 714 701 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 715 702 END_2D … … 728 715 ! 729 716 DO jl = 1, jpl !-- flux in x-direction 730 DO_2D( 1, 0, 1, 0 )717 DO_2D( 1, 1, 1, 0 ) 731 718 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 732 719 END_2D … … 735 722 736 723 DO jl = 1, jpl !-- first guess of tracer from u-flux 737 DO_2D( 0, 0, 0, 0 )724 DO_2D( 1, 1, 0, 0 ) 738 725 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 739 726 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 742 729 END_2D 743 730 END DO 744 #if defined key_mpi3745 CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )746 #else747 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )748 #endif749 731 750 732 DO jl = 1, jpl !-- flux in y-direction 751 DO_2D( 1, 0, 1, 0 )733 DO_2D( 1, 0, 0, 0 ) 752 734 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 753 735 END_2D … … 758 740 ! 759 741 DO jl = 1, jpl !-- flux in y-direction 760 DO_2D( 1, 0, 1, 0)742 DO_2D( 1, 0, 1, 1 ) 761 743 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 762 744 END_2D … … 765 747 ! 766 748 DO jl = 1, jpl !-- first guess of tracer from v-flux 767 DO_2D( 0, 0, 0, 0)749 DO_2D( 0, 0, 1, 1 ) 768 750 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 769 751 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 772 754 END_2D 773 755 END DO 774 #if defined key_mpi3775 CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )776 #else777 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )778 #endif779 756 ! 780 757 DO jl = 1, jpl !-- flux in x-direction 781 DO_2D( 1, 0, 1, 0 )758 DO_2D( 0, 0, 1, 0 ) 782 759 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 783 760 END_2D … … 952 929 ! 953 930 DO jl = 1, jpl 954 DO_2D( 1, 0, 1, 0 )931 DO_2D( 0, 0, 1, 0 ) 955 932 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 956 933 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 961 938 ! 962 939 DO jl = 1, jpl 963 DO_2D( 1, 0, 1, 0 )940 DO_2D( 0, 0, 1, 0 ) 964 941 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 965 942 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 971 948 ! 972 949 DO jl = 1, jpl 973 DO_2D( 1, 0, 1, 0 )950 DO_2D( 0, 0, 1, 0 ) 974 951 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 975 952 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 985 962 ! 986 963 DO jl = 1, jpl 987 DO_2D( 1, 0, 1, 0 )964 DO_2D( 0, 0, 1, 0 ) 988 965 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 989 966 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 999 976 ! 1000 977 DO jl = 1, jpl 1001 DO_2D( 1, 0, 1, 0 )978 DO_2D( 0, 0, 1, 0 ) 1002 979 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 1003 980 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 1020 997 IF( ll_neg ) THEN 1021 998 DO jl = 1, jpl 1022 DO_2D( 1, 0, 1, 0 )999 DO_2D( 0, 0, 1, 0 ) 1023 1000 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1024 1001 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 1030 1007 ! !-- High order flux in i-direction --! 1031 1008 DO jl = 1, jpl 1032 DO_2D( 1, 0, 1, 0 )1009 DO_2D( 0, 0, 1, 0 ) 1033 1010 pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 1034 1011 END_2D … … 1096 1073 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1097 1074 DO jl = 1, jpl 1098 DO_2D( 1, 0, 1, 0 )1075 DO_2D( 1, 0, 0, 0 ) 1099 1076 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1100 1077 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1104 1081 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1105 1082 DO jl = 1, jpl 1106 DO_2D( 1, 0, 1, 0 )1083 DO_2D( 1, 0, 0, 0 ) 1107 1084 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1108 1085 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & … … 1113 1090 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1114 1091 DO jl = 1, jpl 1115 DO_2D( 1, 0, 1, 0 )1092 DO_2D( 1, 0, 0, 0 ) 1116 1093 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1117 1094 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1126 1103 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1127 1104 DO jl = 1, jpl 1128 DO_2D( 1, 0, 1, 0 )1105 DO_2D( 1, 0, 0, 0 ) 1129 1106 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1130 1107 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1139 1116 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1140 1117 DO jl = 1, jpl 1141 DO_2D( 1, 0, 1, 0 )1118 DO_2D( 1, 0, 0, 0 ) 1142 1119 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1143 1120 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1160 1137 IF( ll_neg ) THEN 1161 1138 DO jl = 1, jpl 1162 DO_2D( 1, 0, 1, 0 )1139 DO_2D( 1, 0, 0, 0 ) 1163 1140 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1164 1141 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & … … 1170 1147 ! !-- High order flux in j-direction --! 1171 1148 DO jl = 1, jpl 1172 DO_2D( 1, 0, 1, 0 )1149 DO_2D( 1, 0, 0, 0 ) 1173 1150 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1174 1151 END_2D … … 1206 1183 ! -------------------------------------------------- 1207 1184 DO jl = 1, jpl 1208 DO_2D( 1, 0, 1, 0 )1185 DO_2D( 0, 0, 1, 0 ) 1209 1186 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1187 END_2D 1188 DO_2D( 1, 0, 0, 0 ) 1210 1189 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1211 1190 END_2D … … 1325 1304 ! --------------------------------- 1326 1305 DO jl = 1, jpl 1327 DO_2D( 1, 0, 1, 0 )1306 DO_2D( 0, 0, 1, 0 ) 1328 1307 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1329 1308 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) … … 1336 1315 END_2D 1337 1316 1338 DO_2D( 1, 0, 1, 0 )1317 DO_2D( 1, 0, 0, 0 ) 1339 1318 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1340 1319 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) … … 1709 1688 END SUBROUTINE Hsnow 1710 1689 1690 SUBROUTINE icemax3D( pice , pmax ) 1691 !!--------------------------------------------------------------------- 1692 !! *** ROUTINE icemax3D *** 1693 !! ** Purpose : compute the max of the 9 points around 1694 !!---------------------------------------------------------------------- 1695 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1696 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1697 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1698 INTEGER :: ji, jj, jl ! dummy loop indices 1699 !!---------------------------------------------------------------------- 1700 DO jl = 1, jpl 1701 DO jj = Njs0-1, Nje0+1 1702 DO ji = Nis0, Nie0 1703 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1704 END DO 1705 END DO 1706 DO jj = Njs0, Nje0 1707 DO ji = Nis0, Nie0 1708 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1709 END DO 1710 END DO 1711 END DO 1712 END SUBROUTINE icemax3D 1713 1714 SUBROUTINE icemax4D( pice , pmax ) 1715 !!--------------------------------------------------------------------- 1716 !! *** ROUTINE icemax4D *** 1717 !! ** Purpose : compute the max of the 9 points around 1718 !!---------------------------------------------------------------------- 1719 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1720 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1721 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1722 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1723 !!---------------------------------------------------------------------- 1724 jlay = SIZE( pice , 3 ) ! size of input arrays 1725 DO jl = 1, jpl 1726 DO jk = 1, jlay 1727 DO jj = Njs0-1, Nje0+1 1728 DO ji = Nis0, Nie0 1729 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1730 END DO 1731 END DO 1732 DO jj = Njs0, Nje0 1733 DO ji = Nis0, Nie0 1734 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1735 END DO 1736 END DO 1737 END DO 1738 END DO 1739 END SUBROUTINE icemax4D 1711 1740 1712 1741 #else -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rdgrft.F90
r13630 r13877 349 349 ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 350 350 apartf(ji,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jl-1) ) * & 351 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar 351 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar ) * z1_gstar ) 352 352 ELSE 353 353 apartf(ji,jl) = 0._wp -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rhg_evp.F90
r13807 r13877 148 148 ! 149 149 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 150 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension 150 151 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 151 152 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: … … 170 171 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 171 172 !! --- diags 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 173 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 173 175 !! --- SIMIP diags 174 176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 760 762 & ) * r1_e1e2t(ji,jj) 761 763 zdt2 = zdt * zdt 764 765 zten_i(ji,jj) = zdt 762 766 763 767 ! shear**2 at T points (doc eq. A16) … … 775 779 776 780 ! delta at T points 777 z delta(ji,jj) = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )778 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -z delta(ji,jj)) ) ! 0 if delta=0779 pdelta_i(ji,jj) = z delta(ji,jj) + rn_creepl * rswitch781 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 782 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 783 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 780 784 781 785 END_2D 782 786 #if defined key_mpi3 783 CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 784 #else 785 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp ) 787 CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 788 & zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) 789 #else 790 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, & 791 & zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp ) 786 792 #endif 787 793 788 794 ! --- Store the stress tensor for the next time step --- ! 789 #if defined key_mpi3790 CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )791 #else792 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )793 #endif794 795 pstress1_i (:,:) = zs1 (:,:) 795 796 pstress2_i (:,:) = zs2 (:,:) … … 825 826 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 826 827 827 ! --- stress tensor--- !828 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN829 ! 830 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )828 ! --- Stress tensor invariants (SIMIP diags) --- ! 829 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 830 ! 831 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 831 832 ! 832 DO_2D( 0, 0, 0, 0 ) 833 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 834 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 835 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 836 837 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 838 839 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 840 841 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 842 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 843 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 844 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 845 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 846 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 847 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 848 END_2D 849 #if defined key_mpi3 850 CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 851 #else 852 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp ) 853 #endif 854 ! 855 CALL iom_put( 'isig1' , zsig1 ) 856 CALL iom_put( 'isig2' , zsig2 ) 857 CALL iom_put( 'isig3' , zsig3 ) 858 ! 859 ! Stress tensor invariants (normal and shear stress N/m) 860 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 861 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 862 863 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 833 DO_2D( 1, 1, 1, 1 ) 834 835 ! Ice stresses 836 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 837 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 838 ! I know, this can be confusing... 839 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 840 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 841 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 842 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 843 844 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 845 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 846 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 847 848 END_2D 849 ! 850 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 851 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 852 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 853 854 DEALLOCATE ( zsig_I, zsig_II ) 855 856 ENDIF 857 858 ! --- Normalized stress tensor principal components --- ! 859 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 860 ! Recommendation 1 : we use ice strength, not replacement pressure 861 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 862 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 863 ! 864 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 865 ! 866 DO_2D( 1, 1, 1, 1 ) 867 868 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 869 ! and **deformations** at current iterates 870 ! following Lemieux & Dupont (2020) 871 zfac = zp_delt(ji,jj) 872 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 873 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 874 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 875 876 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 877 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 878 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 879 880 ! Normalized principal stresses (used to display the ellipse) 881 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 882 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 883 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 884 END_2D 885 ! 886 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 887 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 888 889 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 890 864 891 ENDIF 865 892 … … 1008 1035 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 1009 1036 ! close file 1010 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid)1037 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) 1011 1038 ENDIF 1012 1039 -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceitd.F90
r13571 r13877 627 627 END_2D 628 628 ! 629 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 630 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )631 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )632 !633 DO ji = 1, npti634 jdonor(ji,jl) = jl635 ! how much of a_i you send in cat sup is somewhat arbitrary636 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 637 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 638 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 639 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 640 !! zdaice(ji,jl) = a_i_1d(ji) 641 !! zdvice(ji,jl) = v_i_1d(ji)642 !!clem: these are from UCL and work ok 643 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp644 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1)) * 0.5_wp645 END DO646 !647 IF( npti > 0 ) THEN629 IF( npti > 0 ) THEN 630 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 631 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 632 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 633 ! 634 DO ji = 1, npti 635 jdonor(ji,jl) = jl 636 ! how much of a_i you send in cat sup is somewhat arbitrary 637 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 638 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 639 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 640 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 641 !! zdaice(ji,jl) = a_i_1d(ji) 642 !! zdvice(ji,jl) = v_i_1d(ji) 643 !!clem: these are from UCL and work ok 644 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp 645 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 646 END DO 647 ! 648 648 CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl=>jl+1 649 649 ! Reset shift parameters … … 666 666 END_2D 667 667 ! 668 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok669 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok670 DO ji = 1, npti671 jdonor(ji,jl) = jl + 1672 zdaice(ji,jl) = a_i_1d(ji)673 zdvice(ji,jl) = v_i_1d(ji)674 END DO675 !676 668 IF( npti > 0 ) THEN 669 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 670 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 671 DO ji = 1, npti 672 jdonor(ji,jl) = jl + 1 673 zdaice(ji,jl) = a_i_1d(ji) 674 zdvice(ji,jl) = v_i_1d(ji) 675 END DO 676 ! 677 677 CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl+1=>jl 678 678 ! Reset shift parameters -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icestp.F90
r13571 r13877 55 55 USE icedyn ! sea-ice: dynamics 56 56 USE icethd ! sea-ice: thermodynamics 57 USE icecor ! sea-ice: corrections58 57 USE iceupdate ! sea-ice: sea surface boundary condition update 59 58 USE icedia ! sea-ice: budget diagnostics … … 86 85 PUBLIC ice_init ! called by sbcmod.F90 87 86 87 !! * Substitutions 88 # include "do_loop_substitute.h90" 88 89 !!---------------------------------------------------------------------- 89 90 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 160 161 IF( ln_icedyn .AND. .NOT.lk_c1d ) & 161 162 & CALL ice_dyn( kt, Kmm ) ! -- Ice dynamics 163 ! 164 CALL diag_trends( 1 ) ! record dyn trends 162 165 ! 163 166 ! !== lateral boundary conditions ==! … … 188 191 IF( ln_icethd ) CALL ice_thd( kt ) ! -- Ice thermodynamics 189 192 ! 190 CALL ice_cor( kt , 2 ) ! -- Corrections 191 ! 193 CALL diag_trends( 2 ) ! record thermo trends 192 194 CALL ice_var_glo2eqv ! necessary calls (at least for coupling) 193 195 CALL ice_var_agg( 2 ) ! necessary calls (at least for coupling) … … 196 198 ! 197 199 IF( ln_icediahsb ) CALL ice_dia( kt ) ! -- Diagnostics outputs 200 ! 201 IF( ln_icediachk ) CALL ice_drift_wri( kt ) ! -- Diagnostics outputs for conservation 198 202 ! 199 203 CALL ice_wri( kt ) ! -- Ice outputs … … 208 212 ! --- Ocean time step --- ! 209 213 !-------------------------! 210 IF( ln_icedyn ) CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )! -- update surface ocean stresses214 CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) ) ! -- update surface ocean stresses 211 215 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 212 216 ! … … 281 285 ! 282 286 CALL ice_dia_init ! initialization for diags 287 ! 288 CALL ice_drift_init ! initialization for diags of conservation 283 289 ! 284 290 fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction … … 341 347 ENDIF 342 348 ! 343 IF( ln_bdy .AND. ln_icediachk ) CALL ctl_warn('par_init: online conservation check does not work with BDY')344 !345 349 rDt_ice = REAL(nn_fsbc) * rn_Dt !--- sea-ice timestep and its inverse 346 350 r1_Dt_ice = 1._wp / rDt_ice … … 392 396 !! of the time step 393 397 !!---------------------------------------------------------------------- 394 INTEGER :: ji, jj ! dummy loop index 395 !!---------------------------------------------------------------------- 396 sfx (:,:) = 0._wp ; 397 sfx_bri(:,:) = 0._wp ; sfx_lam(:,:) = 0._wp 398 sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp 399 sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp 400 sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp 401 sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp 402 ! 403 wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp 404 wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp 405 wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp 406 wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp 407 wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp 408 wfx_spr(:,:) = 0._wp ; wfx_lam(:,:) = 0._wp 409 wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp 410 wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 411 wfx_snw_sni(:,:) = 0._wp 412 wfx_pnd(:,:) = 0._wp 413 414 hfx_thd(:,:) = 0._wp ; 415 hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp 416 hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp 417 hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp 418 hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp 419 hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp 420 hfx_err_dif(:,:) = 0._wp 421 wfx_err_sub(:,:) = 0._wp 422 ! 423 diag_heat(:,:) = 0._wp ; diag_sice(:,:) = 0._wp 424 diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp 425 426 ! SIMIP diagnostics 427 qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes 428 t_si (:,:,:) = rt0 ! temp at the ice-snow interface 429 430 tau_icebfr (:,:) = 0._wp ! landfast ice param only (clem: important to keep the init here) 431 cnd_ice (:,:,:) = 0._wp ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 432 qcn_ice (:,:,:) = 0._wp ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 433 qtr_ice_bot(:,:,:) = 0._wp ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 434 qsb_ice_bot(:,:) = 0._wp ! (needed if ln_icethd=F) 435 ! 436 ! for control checks (ln_icediachk) 437 diag_trp_vi(:,:) = 0._wp ; diag_trp_vs(:,:) = 0._wp 438 diag_trp_ei(:,:) = 0._wp ; diag_trp_es(:,:) = 0._wp 439 diag_trp_sv(:,:) = 0._wp 440 398 INTEGER :: ji, jj, jl ! dummy loop index 399 !!---------------------------------------------------------------------- 400 401 DO_2D( 1, 1, 1, 1 ) 402 sfx (ji,jj) = 0._wp ; 403 sfx_bri(ji,jj) = 0._wp ; sfx_lam(ji,jj) = 0._wp 404 sfx_sni(ji,jj) = 0._wp ; sfx_opw(ji,jj) = 0._wp 405 sfx_bog(ji,jj) = 0._wp ; sfx_dyn(ji,jj) = 0._wp 406 sfx_bom(ji,jj) = 0._wp ; sfx_sum(ji,jj) = 0._wp 407 sfx_res(ji,jj) = 0._wp ; sfx_sub(ji,jj) = 0._wp 408 ! 409 wfx_snw(ji,jj) = 0._wp ; wfx_ice(ji,jj) = 0._wp 410 wfx_sni(ji,jj) = 0._wp ; wfx_opw(ji,jj) = 0._wp 411 wfx_bog(ji,jj) = 0._wp ; wfx_dyn(ji,jj) = 0._wp 412 wfx_bom(ji,jj) = 0._wp ; wfx_sum(ji,jj) = 0._wp 413 wfx_res(ji,jj) = 0._wp ; wfx_sub(ji,jj) = 0._wp 414 wfx_spr(ji,jj) = 0._wp ; wfx_lam(ji,jj) = 0._wp 415 wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp 416 wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp 417 wfx_snw_sni(ji,jj) = 0._wp 418 wfx_pnd(ji,jj) = 0._wp 419 420 hfx_thd(ji,jj) = 0._wp ; 421 hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp 422 hfx_bog(ji,jj) = 0._wp ; hfx_dyn(ji,jj) = 0._wp 423 hfx_bom(ji,jj) = 0._wp ; hfx_sum(ji,jj) = 0._wp 424 hfx_res(ji,jj) = 0._wp ; hfx_sub(ji,jj) = 0._wp 425 hfx_spr(ji,jj) = 0._wp ; hfx_dif(ji,jj) = 0._wp 426 hfx_err_dif(ji,jj) = 0._wp 427 wfx_err_sub(ji,jj) = 0._wp 428 ! 429 diag_heat(ji,jj) = 0._wp ; diag_sice(ji,jj) = 0._wp 430 diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp 431 432 tau_icebfr (ji,jj) = 0._wp ! landfast ice param only (clem: important to keep the init here) 433 qsb_ice_bot(ji,jj) = 0._wp ! (needed if ln_icethd=F) 434 435 fhld(ji,jj) = 0._wp ! needed if ln_icethd=F 436 437 ! for control checks (ln_icediachk) 438 diag_trp_vi(ji,jj) = 0._wp ; diag_trp_vs(ji,jj) = 0._wp 439 diag_trp_ei(ji,jj) = 0._wp ; diag_trp_es(ji,jj) = 0._wp 440 diag_trp_sv(ji,jj) = 0._wp 441 ! 442 diag_adv_mass(ji,jj) = 0._wp 443 diag_adv_salt(ji,jj) = 0._wp 444 diag_adv_heat(ji,jj) = 0._wp 445 END_2D 446 447 DO jl = 1, jpl 448 DO_2D( 1, 1, 1, 1 ) 449 ! SIMIP diagnostics 450 t_si (ji,jj,jl) = rt0 ! temp at the ice-snow interface 451 qcn_ice_bot(ji,jj,jl) = 0._wp 452 qcn_ice_top(ji,jj,jl) = 0._wp ! conductive fluxes 453 cnd_ice (ji,jj,jl) = 0._wp ! effective conductivity at the top of ice/snow (ln_cndflx=T) 454 qcn_ice (ji,jj,jl) = 0._wp ! conductive flux (ln_cndflx=T & ln_cndemule=T) 455 qtr_ice_bot(ji,jj,jl) = 0._wp ! part of solar radiation transmitted through the ice needed at least for outputs 456 END_2D 457 ENDDO 458 441 459 END SUBROUTINE diag_set0 460 461 462 SUBROUTINE diag_trends( kn ) 463 !!---------------------------------------------------------------------- 464 !! *** ROUTINE diag_trends *** 465 !! 466 !! ** purpose : diagnostics of the trends. Used for conservation purposes 467 !! and outputs 468 !!---------------------------------------------------------------------- 469 INTEGER, INTENT(in) :: kn ! 1 = after dyn ; 2 = after thermo 470 !!---------------------------------------------------------------------- 471 ! 472 ! --- trends of heat, salt, mass (used for conservation controls) 473 IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN 474 ! 475 diag_heat(:,:) = diag_heat(:,:) & 476 & - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice & 477 & - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 478 diag_sice(:,:) = diag_sice(:,:) & 479 & + SUM( sv_i(:,:,:) - sv_i_b(:,:,:) , dim=3 ) * r1_Dt_ice * rhoi 480 diag_vice(:,:) = diag_vice(:,:) & 481 & + SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhoi 482 diag_vsnw(:,:) = diag_vsnw(:,:) & 483 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos 484 ! 485 IF( kn == 2 ) CALL iom_put ( 'hfxdhc' , diag_heat ) ! output of heat trend 486 ! 487 ENDIF 488 ! 489 ! --- trends of concentration (used for simip outputs) 490 IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 491 ! 492 diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 493 ! 494 IF( kn == 1 ) CALL iom_put( 'afxdyn' , diag_aice ) ! dyn trend 495 IF( kn == 2 ) CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) ! thermo trend 496 IF( kn == 2 ) CALL iom_put( 'afxtot' , diag_aice ) ! total trend 497 ! 498 ENDIF 499 ! 500 END SUBROUTINE diag_trends 442 501 443 502 #else -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icetab.F90
r10069 r13877 40 40 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 41 41 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in ) :: tab2d ! input 2D field 42 REAL(wp), DIMENSION(ndim1d,jpl) , INTENT( 42 REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) :: tab1d ! output 1D field 43 43 ! 44 44 INTEGER :: jl, jn, jid, jjd … … 61 61 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 62 62 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: tab2d ! input 2D field 63 REAL(wp), DIMENSION(ndim1d) , INTENT( 63 REAL(wp), DIMENSION(ndim1d) , INTENT(inout) :: tab1d ! output 1D field 64 64 ! 65 65 INTEGER :: jn , jid, jjd … … 80 80 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 81 81 REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab1d ! input 1D field 82 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( 82 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) :: tab2d ! output 2D field 83 83 ! 84 84 INTEGER :: jl, jn, jid, jjd … … 101 101 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 102 102 REAL(wp), DIMENSION(ndim1d) , INTENT(in ) :: tab1d ! input 1D field 103 REAL(wp), DIMENSION(jpi,jpj), INTENT( 103 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: tab2d ! output 2D field 104 104 ! 105 105 INTEGER :: jn , jid, jjd -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd.F90
r13630 r13877 18 18 USE ice ! sea-ice: variables 19 19 !!gm list trop longue ==>>> why not passage en argument d'appel ? 20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot,sprecip, ln_cpl20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 21 21 USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 22 22 & qml_ice, qcn_ice, qtr_ice_top … … 30 30 USE icethd_pnd ! sea-ice: melt ponds 31 31 USE iceitd ! sea-ice: remapping thickness distribution 32 USE icecor ! sea-ice: corrections 32 33 USE icetab ! sea-ice: 1D <==> 2D transformation 33 34 USE icevar ! sea-ice: operations … … 52 53 LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F) 53 54 LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F) 54 LOGICAL :: ln_leadhfx ! 55 LOGICAL :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean 55 56 56 57 !! for convergence tests … … 91 92 ! 92 93 INTEGER :: ji, jj, jk, jl ! dummy loop indices 93 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 94 REAL(wp), PARAMETER :: zfric_umin = 0._wp 95 REAL(wp), PARAMETER :: zch = 0.0057_wp 96 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)94 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 95 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 96 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 97 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 97 98 ! 98 99 !!------------------------------------------------------------------- … … 124 125 & ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 125 126 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 127 zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 128 & ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 126 129 END_2D 127 130 ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean … … 130 133 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 131 134 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 135 zvel(ji,jj) = 0._wp 132 136 END_2D 133 137 ENDIF 134 138 #if defined key_mpi3 135 CALL lbc_lnk_nc_multi( 'icethd', zfric, 'T', 1.0_wp )139 CALL lbc_lnk_nc_multi( 'icethd', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp ) 136 140 #else 137 CALL lbc_lnk ( 'icethd', zfric, 'T',1.0_wp )141 CALL lbc_lnk_multi( 'icethd', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp ) 138 142 #endif 139 143 ! … … 144 148 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 145 149 ! 146 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget147 ! ! practically no "direct lateral ablation"148 !149 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean150 ! ! temperature and turbulent mixing (McPhee, 1992)151 !152 150 ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 153 151 zqld = tmask(ji,jj,1) * rDt_ice * & … … 155 153 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 156 154 157 ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 155 ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 156 ! (mostly<0 but >0 if supercooling) 158 157 zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 159 158 zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0 160 161 ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 159 zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0 160 161 ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 162 ! (mostly>0 but <0 if supercooling) 162 163 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 163 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 164 165 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 164 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 165 166 166 ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 167 167 ! the freezing point, so that we do not have SST < T_freeze 168 ! This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 169 170 !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 171 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 172 173 ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 174 ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 175 IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 176 IF( ln_leadhfx ) THEN ; 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 177 ELSE ; fhld(ji,jj) = 0._wp 178 ENDIF 168 ! This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 169 ! The following formulation is ok for both normal conditions and supercooling 170 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 171 172 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 173 ! qlead is the energy received from the atm. in the leads. 174 ! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2) 175 ! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2) 176 IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 177 ! upper bound for fhld: fhld should be equal to zqld 178 ! but we have to make sure that this heat will not make the sst drop below the freezing point 179 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 180 ! The following formulation is ok for both normal conditions and supercooling 181 fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * 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 182 & - qsb_ice_bot(ji,jj) ) 179 183 qlead(ji,jj) = 0._wp 180 184 ELSE 181 185 fhld (ji,jj) = 0._wp 186 ! upper bound for qlead: qlead should be equal to zqld 187 ! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 188 ! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 189 ! and freezing point is reached if zqfr = zqld - qsb*a/dt 190 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 191 ! The following formulation is ok for both normal conditions and supercooling 192 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 182 193 ENDIF 183 194 ! 184 ! Net heat flux on top of the ice-ocean [W.m-2] 185 ! --------------------------------------------- 186 qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 195 ! If ice is landfast and ice concentration reaches its max 196 ! => stop ice formation in open water 197 IF( zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 ) qlead(ji,jj) = 0._wp 198 ! 199 ! If the grid cell is almost fully covered by ice (no leads) 200 ! => stop ice formation in open water 201 IF( at_i(ji,jj) >= (1._wp - epsi10) ) qlead(ji,jj) = 0._wp 202 ! 203 ! If ln_leadhfx is false 204 ! => do not use energy of the leads to melt sea-ice 205 IF( .NOT.ln_leadhfx ) fhld(ji,jj) = 0._wp 206 ! 187 207 END_2D 188 208 … … 195 215 ENDIF 196 216 197 ! ---------------------------------------------------------------------198 ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]199 ! ---------------------------------------------------------------------200 ! First step here : non solar + precip - qlead - qsensible201 ! Second step in icethd_dh : heat remaining if total melt (zq_rema)202 ! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar203 qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean204 & - qlead(:,:) * r1_Dt_ice & ! heat flux taken from the ocean where there is open water ice formation205 & - at_i (:,:) * qsb_ice_bot(:,:) & ! heat flux taken by sensible flux206 & - at_i (:,:) * fhld (:,:) ! heat flux taken during bottom growth/melt207 ! ! (fhld should be 0 while bott growth)208 217 !-------------------------------------------------------------------------------------------! 209 218 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 259 268 IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- ! 260 269 ! 270 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 271 ! 272 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice ! ice natural aging incrementation 273 ! 261 274 ! convergence tests 262 275 IF( ln_zdf_chkcvg ) THEN … … 422 435 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 423 436 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 424 CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )425 437 ! 426 438 ! ocean surface fields 427 439 CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m ) 428 440 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m ) 441 CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m ) 429 442 ! 430 443 ! to update ice age … … 514 527 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res ) 515 528 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 516 CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d (1:npti), qt_oce_ai )517 529 ! 518 530 CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) ) -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_dh.F90
r13571 r13877 139 139 ! 140 140 DO ji = 1, npti 141 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 141 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) 142 142 zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice ) 143 143 END DO … … 556 556 ! 557 557 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 558 qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice558 !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 559 559 560 560 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_do.F90
r13630 r13877 131 131 132 132 ! Default new ice thickness 133 WHERE( qlead(:,:) < 0._wp .AND. tau_icebfr(:,:) == 0._wp ) ; ht_i_new(:,:) = rn_hinew ! if cooling and no landfast 134 ELSEWHERE ; ht_i_new(:,:) = 0._wp 133 WHERE( qlead(:,:) < 0._wp ) ! cooling 134 ht_i_new(:,:) = rn_hinew 135 ELSEWHERE 136 ht_i_new(:,:) = 0._wp 135 137 END WHERE 136 138 … … 146 148 ! 147 149 DO_2D( 0, 0, 0, 0 ) 148 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast150 IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 149 151 ! -- Wind stress -- ! 150 152 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & … … 202 204 ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 203 205 !------------------------------------------------------------------------------! 204 ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice206 ! it occurs if cooling 205 207 206 208 ! Identify grid points where new ice forms 207 209 npti = 0 ; nptidx(:) = 0 208 210 DO_2D( 1, 1, 1, 1 ) 209 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp) THEN211 IF ( qlead(ji,jj) < 0._wp ) THEN 210 212 npti = npti + 1 211 213 nptidx( npti ) = (jj - 1) * jpi + ji -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceupdate.F90
r13630 r13877 24 24 USE traqsr ! add penetration of solar flux in the calculation of heat budget 25 25 USE icectl ! sea-ice: control prints 26 USE bdy_oce , ONLY : ln_bdy27 26 USE zdfdrg , ONLY : ln_drgice_imp 28 27 ! … … 92 91 ! 93 92 INTEGER :: ji, jj, jl, jk ! dummy loop indices 94 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2)95 93 REAL(wp) :: zqsr ! New solar flux received by the ocean 96 94 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 103 101 WRITE(numout,*)'~~~~~~~~~~~~~~' 104 102 ENDIF 103 104 ! Net heat flux on top of the ice-ocean (W.m-2) 105 !---------------------------------------------- 106 qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 105 107 106 108 ! --- case we bypass ice thermodynamics --- ! … … 115 117 DO_2D( 1, 1, 1, 1 ) 116 118 117 ! Solar heat flux reaching the ocean = zqsr (W.m-2)119 ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2) 118 120 !--------------------------------------------------- 119 121 zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) ) … … 121 123 ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2) 122 124 !--------------------------------------------------- 123 zqmass = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 124 qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr 125 126 ! Add the residual from heat diffusion equation and sublimation (W.m-2) 127 !---------------------------------------------------------------------- 128 qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) + & 129 & ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) 130 125 qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) & 126 & - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) & 127 & + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) & 128 & + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj) 129 131 130 ! New qsr and qns used to compute the oceanic heat flux at the next time step 132 131 !---------------------------------------------------------------------------- 133 qsr(ji,jj) = zqsr 132 ! if warming and some ice remains, then we suppose that the whole solar flux has been consumed to melt the ice 133 ! else ( cooling or no ice left ), then we suppose that no solar flux has been consumed 134 ! 135 IF( fhld(ji,jj) > 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN !-- warming and some ice remains 136 ! solar flux transmitted thru the 1st level of the ocean (i.e. not used by sea-ice) 137 qsr(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * ( 1._wp - frq_m(ji,jj) ) & 138 ! + solar flux transmitted thru ice and the 1st ocean level (also not used by sea-ice) 139 & + SUM( a_i_b(ji,jj,:) * qtr_ice_bot(ji,jj,:) ) * ( 1._wp - frq_m(ji,jj) ) 140 ! 141 ELSE !-- cooling or no ice left 142 qsr(ji,jj) = zqsr 143 ENDIF 144 ! 145 ! the non-solar is simply derived from the solar flux 134 146 qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr 135 147 136 148 ! Mass flux at the atm. surface 137 149 !----------------------------------- … … 140 152 ! Mass flux at the ocean surface 141 153 !------------------------------------ 142 ! case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED) 143 ! ------------------------------------------------------------------------------------- 144 ! The idea of this approach is that the system that we consider is the ICE-OCEAN system 145 ! Thus FW flux = External ( E-P+snow melt) 146 ! Salt flux = Exchanges in the ice-ocean system then converted into FW 147 ! Associated to Ice formation AND Ice melting 148 ! Even if i see Ice melting as a FW and SALT flux 149 ! 150 ! mass flux from ice/ocean 154 ! ice-ocean mass flux 151 155 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 152 156 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 153 154 ! add the snow melt water to snow mass flux to the ocean157 158 ! snw-ocean mass flux 155 159 wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj) 156 157 ! mass flux at the ocean/ice interface 158 fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) ) ! F/M mass flux save at least for biogeochemical model 159 emp(ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) 160 160 161 ! total mass flux at the ocean/ice interface 162 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model 163 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux 161 164 162 165 ! Salt flux at the ocean surface … … 262 265 CALL iom_put ('hfxdif' , hfx_dif ) ! heat flux used for ice temperature change 263 266 CALL iom_put ('hfxsnw' , hfx_snw ) ! heat flux used for snow melt 264 CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion (included in qt_oce_ai)267 CALL iom_put ('hfxerr' , hfx_err_dif ) ! heat flux error after heat diffusion 265 268 266 269 ! heat fluxes associated with mass exchange (freeze/melt/precip...) … … 279 282 !--------- 280 283 #if ! defined key_agrif 281 IF( ln_icediachk .AND. .NOT. ln_bdy) CALL ice_cons_final('iceupdate') ! conservation284 IF( ln_icediachk ) CALL ice_cons_final('iceupdate') ! conservation 282 285 #endif 283 IF( ln_icectl 284 IF( sn_cfctl%l_prtctl 285 IF( ln_timing 286 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 287 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('iceupdate') ! prints 288 IF( ln_timing ) CALL timing_stop ('ice_update') ! timing 286 289 ! 287 290 END SUBROUTINE ice_update_flx -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/BDY/bdyice.F90
r13630 r13877 61 61 !!---------------------------------------------------------------------- 62 62 ! controls 63 IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing 64 IF( ln_icediachk ) CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 65 IF( ln_icediachk ) CALL ice_cons2D (0,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 63 IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing 66 64 ! 67 65 CALL ice_var_glo2eqv … … 120 118 ! 121 119 ! controls 122 IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints 123 IF( ln_icediachk ) CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 124 IF( ln_icediachk ) CALL ice_cons2D (1,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 125 IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing 120 IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints 121 IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing 126 122 ! 127 123 END SUBROUTINE bdy_ice -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DIU/diu_coolskin.F90
r13295 r13877 95 95 !!---------------------------------------------------------------------- 96 96 ! 97 IF( .NOT. ln_blk) CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing")97 IF( .NOT. (ln_blk .OR. ln_abl) ) CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing") 98 98 ! 99 99 DO_2D( 1, 1, 1, 1 ) -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/IOM/iom.F90
r13807 r13877 1935 1935 IF( iom_use(cdname) ) THEN 1936 1936 #if defined key_iomput 1937 IF( SIZE(pfield2d, dim=1) == jpi .AND. SIZE(pfield2d, dim=2) == jpj ) THEN 1938 CALL xios_send_field( cdname, pfield2d(Nis0:Nie0, Njs0:Nje0) ) ! this extraction will create a copy of pfield2d 1939 ELSE 1940 CALL xios_send_field( cdname, pfield2d ) 1941 ENDIF 1937 CALL xios_send_field( cdname, pfield2d ) 1942 1938 #else 1943 1939 WRITE(numout,*) pfield2d ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings … … 1951 1947 IF( iom_use(cdname) ) THEN 1952 1948 #if defined key_iomput 1953 IF( SIZE(pfield2d, dim=1) == jpi .AND. SIZE(pfield2d, dim=2) == jpj ) THEN 1954 CALL xios_send_field( cdname, pfield2d(Nis0:Nie0, Njs0:Nje0) ) ! this extraction will create a copy of pfield2d 1955 ELSE 1956 CALL xios_send_field( cdname, pfield2d ) 1957 ENDIF 1949 CALL xios_send_field( cdname, pfield2d ) 1958 1950 #else 1959 1951 WRITE(numout,*) pfield2d ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings … … 1967 1959 IF( iom_use(cdname) ) THEN 1968 1960 #if defined key_iomput 1969 IF( SIZE(pfield3d, dim=1) == jpi .AND. SIZE(pfield3d, dim=2) == jpj ) THEN 1970 CALL xios_send_field( cdname, pfield3d(Nis0:Nie0, Njs0:Nje0,:) ) ! this extraction will create a copy of pfield3d 1971 ELSE 1972 CALL xios_send_field( cdname, pfield3d ) 1973 ENDIF 1961 CALL xios_send_field( cdname, pfield3d ) 1974 1962 #else 1975 1963 WRITE(numout,*) pfield3d ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings … … 1983 1971 IF( iom_use(cdname) ) THEN 1984 1972 #if defined key_iomput 1985 IF( SIZE(pfield3d, dim=1) == jpi .AND. SIZE(pfield3d, dim=2) == jpj ) THEN 1986 CALL xios_send_field( cdname, pfield3d(Nis0:Nie0, Njs0:Nje0,:) ) ! this extraction will create a copy of pfield3d 1987 ELSE 1988 CALL xios_send_field( cdname, pfield3d ) 1989 ENDIF 1973 CALL xios_send_field( cdname, pfield3d ) 1990 1974 #else 1991 1975 WRITE(numout,*) pfield3d ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings … … 1999 1983 IF( iom_use(cdname) ) THEN 2000 1984 #if defined key_iomput 2001 IF( SIZE(pfield4d, dim=1) == jpi .AND. SIZE(pfield4d, dim=2) == jpj ) THEN 2002 CALL xios_send_field( cdname, pfield4d(Nis0:Nie0, Njs0:Nje0,:,:) ) ! this extraction will create a copy of pfield4d 2003 ELSE 2004 CALL xios_send_field (cdname, pfield4d ) 2005 ENDIF 1985 CALL xios_send_field (cdname, pfield4d ) 2006 1986 #else 2007 1987 WRITE(numout,*) pfield4d ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings … … 2015 1995 IF( iom_use(cdname) ) THEN 2016 1996 #if defined key_iomput 2017 IF( SIZE(pfield4d, dim=1) == jpi .AND. SIZE(pfield4d, dim=2) == jpj ) THEN 2018 CALL xios_send_field( cdname, pfield4d(Nis0:Nie0, Njs0:Nje0,:,:) ) ! this extraction will create a copy of pfield4d 2019 ELSE 2020 CALL xios_send_field (cdname, pfield4d ) 2021 ENDIF 1997 CALL xios_send_field (cdname, pfield4d ) 2022 1998 #else 2023 1999 WRITE(numout,*) pfield4d ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings … … 2225 2201 ! 2226 2202 CALL iom_set_domain_attr("grid_"//cdgrd, ni_glo=Ni0glo,nj_glo=Nj0glo,ibegin=mig0(Nis0)-1,jbegin=mjg0(Njs0)-1,ni=Ni_0,nj=Nj_0) 2227 CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = 0, data_ni = Ni_0, data_jbegin = 0, data_nj = Nj_0)2203 CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj) 2228 2204 !don't define lon and lat for restart reading context. 2229 2205 IF ( .NOT.ldrxios ) & … … 2328 2304 CALL dom_ngb( 180.0_wp, 90.0_wp, ix, iy, 'T' ) ! i-line that passes near the North Pole : Reference latitude (used in plots) 2329 2305 CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig0(Nis0)-1, jbegin=mjg0(Njs0)-1, ni=Ni_0, nj=Nj_0) 2330 CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = 0, data_ni = Ni_0, data_jbegin = 0, data_nj = Nj_0)2306 CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj) 2331 2307 CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp), & 2332 2308 & latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp)) 2333 CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj _0)2309 CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo) 2334 2310 ! 2335 2311 CALL iom_update_file_name('ptr') -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/LBC/lib_mpp.F90
r13571 r13877 516 516 ALLOCATE(todelay(idvar)%y1d(isz)) 517 517 todelay(idvar)%y1d(:) = CMPLX(todelay(idvar)%z1d(:), 0., wp) ! create %y1d, complex variable needed by mpi_sumdd 518 ndelayid(idvar) = MPI_REQUEST_NULL ! initialised request to a valid value 518 519 END IF 519 520 ENDIF … … 523 524 ALLOCATE(todelay(idvar)%z1d(isz), todelay(idvar)%y1d(isz)) ! allocate also %z1d as used for the restart 524 525 CALL mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr ) ! get %y1d 525 todelay(idvar)%z1d(:) = REAL(todelay(idvar)%y1d(:), wp) ! define %z1d from %y1d526 ENDIF 527 528 IF( ndelayid(idvar) > 0 )CALL mpp_delay_rcv( idvar ) ! make sure %z1d is received526 ndelayid(idvar) = MPI_REQUEST_NULL 527 ENDIF 528 529 CALL mpp_delay_rcv( idvar ) ! make sure %z1d is received 529 530 530 531 ! send back pout from todelay(idvar)%z1d defined at previous call … … 535 536 IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 536 537 CALL mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr ) 537 ndelayid(idvar) = 1538 ndelayid(idvar) = MPI_REQUEST_NULL 538 539 IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 539 540 # else … … 596 597 DEALLOCATE(todelay(idvar)%z1d) 597 598 ndelayid(idvar) = -1 ! do as if we had no restart 599 ELSE 600 ndelayid(idvar) = MPI_REQUEST_NULL 598 601 END IF 599 602 ENDIF … … 603 606 ALLOCATE(todelay(idvar)%z1d(isz)) 604 607 CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_DOUBLE_PRECISION, mpi_max, ilocalcomm, ierr ) ! get %z1d 605 ENDIF 606 607 IF( ndelayid(idvar) > 0 ) CALL mpp_delay_rcv( idvar ) ! make sure %z1d is received 608 ndelayid(idvar) = MPI_REQUEST_NULL 609 ENDIF 610 611 CALL mpp_delay_rcv( idvar ) ! make sure %z1d is received 608 612 609 613 ! send back pout from todelay(idvar)%z1d defined at previous call … … 611 615 612 616 ! send p_in into todelay(idvar)%z1d with a non-blocking communication 617 ! (PM) Should we get rid of MPI2 option ? MPI3 was release in 2013. Who is still using MPI2 ? 613 618 # if defined key_mpi2 614 619 IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 615 CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ndelayid(idvar),ierr )620 CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ierr ) 616 621 IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 617 622 # else … … 636 641 !!---------------------------------------------------------------------- 637 642 #if defined key_mpp_mpi 638 IF( ndelayid(kid) /= -2 ) THEN 639 #if ! defined key_mpi2 640 IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 641 CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr ) ! make sure todelay(kid) is received 642 IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.) 643 #endif 644 IF( ASSOCIATED(todelay(kid)%y1d) ) todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp) ! define %z1d from %y1d 645 ndelayid(kid) = -2 ! add flag to know that mpi_wait was already called on kid 646 ENDIF 643 IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.) 644 ! test on ndelayid(kid) useless as mpi_wait return immediatly if the request handle is MPI_REQUEST_NULL 645 CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr ) ! after this ndelayid(kid) = MPI_REQUEST_NULL 646 IF( ln_timing ) CALL tic_tac( .FALSE., ld_global = .TRUE.) 647 IF( ASSOCIATED(todelay(kid)%y1d) ) todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp) ! define %z1d from %y1d 647 648 #endif 648 649 END SUBROUTINE mpp_delay_rcv -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcfwb.F90
r13630 r13877 94 94 snwice_mass_b(:,:) = 0.e0 ! no sea-ice model is being used : no snow+ice mass 95 95 snwice_mass (:,:) = 0.e0 96 snwice_fmass (:,:) = 0.e0 96 97 #endif 97 98 ! -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcmod.F90
r13630 r13877 229 229 CASE DEFAULT !- not supported 230 230 END SELECT 231 IF( ln_diurnal .AND. .NOT. ln_blk) CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" )231 IF( ln_diurnal .AND. .NOT. (ln_blk.OR.ln_abl) ) CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" ) 232 232 ! 233 233 ! !** allocate and set required variables -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/stpctl.F90
r13571 r13877 67 67 REAL(wp) :: zzz ! local real 68 68 REAL(wp), DIMENSION(9) :: zmax, zmaxlocal 69 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 69 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 70 70 LOGICAL, DIMENSION(jpi,jpj,jpk) :: llmsk 71 71 CHARACTER(len=20) :: clname … … 125 125 ! 126 126 llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp ! define only the inner domain 127 ! 128 ll_0oce = .NOT. ANY( llmsk(:,:,1) ) ! no ocean point in the inner domain? 129 ! 127 130 IF( ll_wd ) THEN 128 131 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) ) ! ssh max … … 149 152 ENDIF 150 153 zmax(9) = REAL( nstop, wp ) ! stop indicator 154 ! 151 155 ! !== get global extrema ==! 152 156 ! !== done by all processes if writting run.stat ==! 153 157 IF( ll_colruns ) THEN 154 158 zmaxlocal(:) = zmax(:) 155 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 159 CALL mpp_max( "stpctl", zmax ) ! max over the global domain: ok even of ll_0oce = .true. 156 160 nstop = NINT( zmax(9) ) ! update nstop indicator (now sheared among all local domains) 157 ENDIF 161 ELSE 162 ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 163 IF( ll_0oce ) zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /) ! default "valid" values... 164 ENDIF 165 ! 166 zmax(3) = -zmax(3) ! move back from max(-zz) to min(zz) : easier to manage! 167 zmax(5) = -zmax(5) ! move back from max(-zz) to min(zz) : easier to manage! 168 IF( ll_colruns ) THEN 169 zmaxlocal(3) = -zmaxlocal(3) ! move back from max(-zz) to min(zz) : easier to manage! 170 zmaxlocal(5) = -zmaxlocal(5) ! move back from max(-zz) to min(zz) : easier to manage! 171 ENDIF 172 ! 158 173 ! !== write "run.stat" files ==! 159 174 ! !== done only by 1st subdomain at writting timestep ==! 160 175 IF( ll_wrtruns ) THEN 161 WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 162 istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 163 istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 164 istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 165 istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 166 istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 167 istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 168 IF( ln_zad_Aimp ) THEN 169 istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 170 istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 171 ENDIF 176 WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4) 177 DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) ) 178 istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 179 END DO 172 180 IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 173 181 END IF … … 175 183 ! !== done by all processes at every time step ==! 176 184 ! 177 IF( 178 & 179 & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity180 & 181 & 182 & 183 & 185 IF( zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 186 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 187 & zmax(3) <= 0._wp .OR. & ! negative or zero sea surface salinity 188 & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) 189 & zmax(4) < 0._wp .OR. & ! too large sea surface salinity (keep this line for sea-ice) 190 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR. & ! NaN encounter in the tests 191 & ABS( zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN ! Infinity encounter in the tests 184 192 ! 185 193 iloc(:,:) = 0 … … 221 229 ! 222 230 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 223 CALL wrt_line( ctmp2, kt, '|ssh| max', 224 CALL wrt_line( ctmp3, kt, '|U| max', 225 CALL wrt_line( ctmp4, kt, 'Sal min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )226 CALL wrt_line( ctmp5, kt, 'Sal max', 231 CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 232 CALL wrt_line( ctmp3, kt, '|U| max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 233 CALL wrt_line( ctmp4, kt, 'Sal min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 234 CALL wrt_line( ctmp5, kt, 'Sal max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 227 235 IF( Agrif_Root() ) THEN 228 236 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort* files' -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/SAS/stpctl.F90
r13571 r13877 68 68 REAL(wp) :: zzz ! local real 69 69 REAL(wp), DIMENSION(4) :: zmax, zmaxlocal 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 71 71 LOGICAL, DIMENSION(jpi,jpj) :: llmsk 72 72 CHARACTER(len=20) :: clname … … 124 124 llmsk(:,Nje1: jpj) = .FALSE. 125 125 ! 126 llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp ! test only the inner domain 127 IF( COUNT( llmsk(:,:) ) > 0 ) THEN ! avoid huge values sent back for land processors... 128 zmax(1) = MAXVAL( vt_i (:,:) , mask = llmsk ) ! max ice thickness 129 zmax(2) = MAXVAL( ABS( u_ice(:,:) ) , mask = llmsk ) ! max ice velocity (zonal only) 130 zmax(3) = MAXVAL( -tm_i (:,:) + rt0, mask = llmsk ) ! min ice temperature (in degC) 131 ELSE 132 IF( ll_colruns ) THEN ! default value: must not be kept when calling mpp_max -> must be as small as possible 133 zmax(1:3) = -HUGE(1._wp) 134 ELSE ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 135 zmax(1:3) = 0._wp 136 ENDIF 137 ENDIF 138 zmax(4) = REAL( nstop, wp ) ! stop indicator 126 llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp ! test only the inner domain 127 ! 128 ll_0oce = .NOT. ANY( llmsk(:,:) ) ! no ocean point in the inner domain? 129 ! 130 zmax(1) = MAXVAL( vt_i (:,:) , mask = llmsk ) ! max ice thickness 131 zmax(2) = MAXVAL( ABS( u_ice(:,:) ) , mask = llmsk ) ! max ice velocity (zonal only) 132 zmax(3) = MAXVAL( -tm_i (:,:) + rt0, mask = llmsk ) ! min ice temperature (in degC) 133 zmax(4) = REAL( nstop, wp ) ! stop indicator 134 ! 139 135 ! !== get global extrema ==! 140 136 ! !== done by all processes if writting run.stat ==! 141 137 IF( ll_colruns ) THEN 142 138 zmaxlocal(:) = zmax(:) 143 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 139 CALL mpp_max( "stpctl", zmax ) ! max over the global domain: ok even of ll_0oce = .true. 144 140 nstop = NINT( zmax(4) ) ! update nstop indicator (now sheared among all local domains) 145 ENDIF 141 ELSE 142 ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 143 IF( ll_0oce ) zmax(1:3) = 0._wp ! default "valid" values... 144 ENDIF 145 ! 146 zmax(3) = -zmax(3) ! move back from max(-zz) to min(zz) : easier to manage! 147 IF( ll_colruns ) zmaxlocal(3) = -zmaxlocal(3) ! move back from max(-zz) to min(zz) : easier to manage! 148 ! 146 149 ! !== write "run.stat" files ==! 147 150 ! !== done only by 1st subdomain at writting timestep ==! 148 151 IF( ll_wrtruns ) THEN 149 WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3)150 istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )151 istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )152 istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) )152 WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3) 153 DO ji = 1, 3 154 istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 155 END DO 153 156 IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 154 157 END IF … … 158 161 IF( zmax(1) > 100._wp .OR. & ! too large ice thickness maximum ( > 100 m) 159 162 & zmax(2) > 10._wp .OR. & ! too large ice velocity ( > 10 m/s) 160 & zmax(3) >101._wp .OR. & ! too cold ice temperature ( < -100 degC)163 & zmax(3) < -101._wp .OR. & ! too cold ice temperature ( < -100 degC) 161 164 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR. & ! NaN encounter in the tests 162 165 & ABS( zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN ! Infinity encounter in the tests … … 192 195 ! 193 196 WRITE(ctmp1,*) ' stp_ctl: ice_thick > 100 m or |ice_vel| > 10 m/s or ice_temp < -100 degC or NaN encounter in the tests' 194 CALL wrt_line( ctmp2, kt, 'ice_thick max', 195 CALL wrt_line( ctmp3, kt, '|ice_vel| max', 196 CALL wrt_line( ctmp4, kt, 'ice_temp min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )197 CALL wrt_line( ctmp2, kt, 'ice_thick max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 198 CALL wrt_line( ctmp3, kt, '|ice_vel| max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 199 CALL wrt_line( ctmp4, kt, 'ice_temp min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 197 200 IF( Agrif_Root() ) THEN 198 201 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort* files' -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CANAL/MY_SRC/stpctl.F90
r13632 r13877 67 67 REAL(wp) :: zzz ! local real 68 68 REAL(wp), DIMENSION(9) :: zmax, zmaxlocal 69 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 69 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 70 70 LOGICAL, DIMENSION(jpi,jpj,jpk) :: llmsk 71 71 CHARACTER(len=20) :: clname … … 125 125 ! 126 126 llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp ! define only the inner domain 127 ! 128 ll_0oce = .NOT. ANY( llmsk(:,:,1) ) ! no ocean point in the inner domain? 129 ! 127 130 IF( ll_wd ) THEN 128 131 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) ) ! ssh max … … 149 152 ENDIF 150 153 zmax(9) = REAL( nstop, wp ) ! stop indicator 154 ! 151 155 ! !== get global extrema ==! 152 156 ! !== done by all processes if writting run.stat ==! 153 157 IF( ll_colruns ) THEN 154 158 zmaxlocal(:) = zmax(:) 155 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 159 CALL mpp_max( "stpctl", zmax ) ! max over the global domain: ok even of ll_0oce = .true. 156 160 nstop = NINT( zmax(9) ) ! update nstop indicator (now sheared among all local domains) 157 ENDIF 161 ELSE 162 ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 163 IF( ll_0oce ) zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /) ! default "valid" values... 164 ENDIF 165 ! 166 zmax(3) = -zmax(3) ! move back from max(-zz) to min(zz) : easier to manage! 167 zmax(5) = -zmax(5) ! move back from max(-zz) to min(zz) : easier to manage! 168 IF( ll_colruns ) THEN 169 zmaxlocal(3) = -zmaxlocal(3) ! move back from max(-zz) to min(zz) : easier to manage! 170 zmaxlocal(5) = -zmaxlocal(5) ! move back from max(-zz) to min(zz) : easier to manage! 171 ENDIF 172 ! 158 173 ! !== write "run.stat" files ==! 159 174 ! !== done only by 1st subdomain at writting timestep ==! 160 175 IF( ll_wrtruns ) THEN 161 WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 162 istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 163 istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 164 istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 165 istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 166 istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 167 istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 168 IF( ln_zad_Aimp ) THEN 169 istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 170 istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 171 ENDIF 176 WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4) 177 DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) ) 178 istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 179 END DO 172 180 IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 173 181 END IF … … 175 183 ! !== done by all processes at every time step ==! 176 184 ! 177 IF( 178 & 179 !!$ & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity180 !!$ & 181 !!$ & 182 & 183 & 185 IF( zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 186 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 187 !!$ & zmax(3) <= 0._wp .OR. & ! negative or zero sea surface salinity 188 !!$ & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) 189 !!$ & zmax(4) < 0._wp .OR. & ! too large sea surface salinity (keep this line for sea-ice) 190 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR. & ! NaN encounter in the tests 191 & ABS( zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN ! Infinity encounter in the tests 184 192 ! 185 193 iloc(:,:) = 0 … … 207 215 ELSE ! find local min and max locations: 208 216 ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc 209 llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp 217 llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp ! define only the inner domain 210 218 iloc(1:2,1) = MAXLOC( ABS( ssh(:,:, Kmm)), mask = llmsk(:,:,1) ) 211 219 llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp ! define only the inner domain … … 221 229 ! 222 230 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 223 CALL wrt_line( ctmp2, kt, '|ssh| max', 224 CALL wrt_line( ctmp3, kt, '|U| max', 225 CALL wrt_line( ctmp4, kt, 'Sal min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )226 CALL wrt_line( ctmp5, kt, 'Sal max', 231 CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 232 CALL wrt_line( ctmp3, kt, '|U| max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 233 CALL wrt_line( ctmp4, kt, 'Sal min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 234 CALL wrt_line( ctmp5, kt, 'Sal max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 227 235 IF( Agrif_Root() ) THEN 228 236 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort* files' -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CPL_OASIS/EXPREF/file_def_nemo-ice.xml
r12663 r13877 53 53 <field field_ref="normstr" name="normstr" /> 54 54 <field field_ref="sheastr" name="sheastr" /> 55 <field field_ref="isig1" name="isig1" /> 56 <field field_ref="isig2" name="isig2" /> 57 <field field_ref="isig3" name="isig3" /> 55 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 56 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 58 57 59 58 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_ADV2D/EXPREF/file_def_nemo-ice.xml
r10516 r13877 55 55 <field field_ref="normstr" name="normstr" /> 56 56 <field field_ref="sheastr" name="sheastr" /> 57 <field field_ref="isig1" name="isig1" /> 58 <field field_ref="isig2" name="isig2" /> 59 <field field_ref="isig3" name="isig3" /> 57 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 58 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 60 59 61 60 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_AGRIF/EXPREF/file_def_nemo-ice.xml
r11159 r13877 53 53 <field field_ref="normstr" name="normstr" /> 54 54 <field field_ref="sheastr" name="sheastr" /> 55 <field field_ref="isig1" name="isig1" /> 56 <field field_ref="isig2" name="isig2" /> 57 <field field_ref="isig3" name="isig3" /> 55 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 56 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 58 57 59 58 <!-- heat fluxes --> -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/EXPREF/namelist_cfg
r13632 r13877 477 477 / 478 478 !----------------------------------------------------------------------- 479 !-----------------------------------------------------------------------480 /481 !-----------------------------------------------------------------------482 479 &namhsb ! Heat and salt budgets (default: OFF) 483 480 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/dtatsd.F90
r13295 r13877 191 191 ENDIF 192 192 ! 193 DO_2D( 1, 1, 1, 1 ) 193 DO_2D( 1, 1, 1, 1 ) ! vertical interpolation of T & S 194 194 DO jk = 1, jpk ! determines the intepolated T-S profiles at each (i,j) points 195 195 zl = gdept_0(ji,jj,jk) -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/eosbn2.F90
r13295 r13877 182 182 !! * Substitutions 183 183 # include "do_loop_substitute.h90" 184 # include "domzgr_substitute.h90" 184 185 !!---------------------------------------------------------------------- 185 186 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 965 966 IF( ln_timing ) CALL timing_start('bn2') 966 967 ! 967 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 968 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90 968 969 zrw = ( gdepw(ji,jj,jk ,Kmm) - gdept(ji,jj,jk,Kmm) ) & 969 970 & / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) … … 1723 1724 ! 1724 1725 CASE( np_leos ) !== Linear ISOMIP EOS ==! 1726 1727 r1_S0 = 0.875_wp/35.16504_wp ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct) 1728 1725 1729 IF(lwp) THEN 1726 1730 WRITE(numout,*) … … 1731 1735 WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0 1732 1736 ENDIF 1737 l_useCT = .TRUE. ! Use conservative temperature 1733 1738 ! 1734 1739 CASE DEFAULT !== ERROR in neos ==! -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isf_oce.F90
r12077 r13877 75 75 ! 76 76 ! 2.1 -------- ice shelf cavity parameter -------------- 77 LOGICAL , PUBLIC :: l_isfoasis 77 LOGICAL , PUBLIC :: l_isfoasis = .FALSE. 78 78 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfload !: ice shelf load 79 79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_oasis -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfcavgam.F90
r12905 r13877 30 30 PUBLIC isfcav_gammats 31 31 32 # include "domzgr_substitute.h90" 32 33 !!---------------------------------------------------------------------- 33 34 !! NEMO/OCE 4.0 , NEMO Consortium (2018) -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfstp.F90
r12905 r13877 13 13 !! isfstp : compute iceshelf melt and heat flux 14 14 !!---------------------------------------------------------------------- 15 !16 15 USE isf_oce ! isf variables 17 16 USE isfload, ONLY: isf_load ! ice shelf load … … 21 20 USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables 22 21 23 USE dom_oce, ONLY: ht, e3t, ln_isfcav, ln_linssh ! ocean space and time domain 22 USE dom_oce ! ocean space and time domain 23 USE oce , ONLY: ssh ! sea surface height 24 24 USE domvvl, ONLY: ln_vvl_zstar ! zstar logical 25 25 USE zdfdrg, ONLY: r_Cdmin_top, r_ke0_top ! vertical physics: top/bottom drag coef. … … 31 31 32 32 IMPLICIT NONE 33 34 33 PRIVATE 35 34 36 35 PUBLIC isf_stp, isf_init, isf_nam ! routine called in sbcmod and divhor 37 36 37 !! * Substitutions 38 # include "domzgr_substitute.h90" 38 39 !!---------------------------------------------------------------------- 39 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 41 42 !! Software governed by the CeCILL license (see ./LICENSE) 42 43 !!---------------------------------------------------------------------- 44 43 45 CONTAINS 44 46 … … 60 62 INTEGER, INTENT(in) :: kt ! ocean time step 61 63 INTEGER, INTENT(in) :: Kmm ! ocean time level index 64 !!---------------------------------------------------------------------- 65 INTEGER :: jk ! loop index 66 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t ! e3t 62 67 !!--------------------------------------------------------------------- 63 68 ! … … 78 83 ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 79 84 rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:) 80 CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 85 DO jk = 1, jpk 86 ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 87 END DO 88 CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 81 89 ! 82 90 ! 1.3: compute ice shelf melt … … 100 108 ! by simplicity, we assume the top level where param applied do not change with time (done in init part) 101 109 rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:) 102 CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 110 DO jk = 1, jpk 111 ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 112 END DO 113 CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 103 114 ! 104 115 ! 2.3: compute ice shelf melt -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/istate.F90
r13295 r13877 24 24 USE dom_oce ! ocean space and time domain 25 25 USE daymod ! calendar 26 USE divhor ! horizontal divergence (div_hor routine)27 26 USE dtatsd ! data temperature and salinity (dta_tsd routine) 28 27 USE dtauvd ! data: U & V current (dta_uvd routine) … … 35 34 USE lib_mpp ! MPP library 36 35 USE restart ! restart 36 #if defined key_agrif 37 USE agrif_oce_interp 38 USE agrif_oce 39 #endif 37 40 38 41 IMPLICIT NONE … … 43 46 !! * Substitutions 44 47 # include "do_loop_substitute.h90" 48 # include "domzgr_substitute.h90" 45 49 !!---------------------------------------------------------------------- 46 50 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 59 63 ! 60 64 INTEGER :: ji, jj, jk ! dummy loop indices 65 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgdept ! 3D table !!st patch to use gdept subtitute 61 66 !!gm see comment further down 62 67 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: zuvd ! U & V data workspace … … 70 75 !!gm Why not include in the first call of dta_tsd ? 71 76 !!gm probably associated with the use of internal damping... 72 77 CALL dta_tsd_init ! Initialisation of T & S input data 73 78 !!gm to be moved in usrdef of C1D case 74 79 ! IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data … … 84 89 #endif 85 90 91 #if defined key_agrif 92 IF ( (.NOT.Agrif_root()).AND.ln_init_chfrpar ) THEN 93 numror = 0 ! define numror = 0 -> no restart file to read 94 ln_1st_euler = .true. ! Set time-step indicator at nit000 (euler forward) 95 CALL day_init 96 CALL agrif_istate( Kbb, Kmm, Kaa ) ! Interp from parent 97 ! 98 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) 99 ssh (:,:,Kmm) = ssh(:,:,Kbb) 100 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 101 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 102 ELSE 103 #endif 86 104 IF( ln_rstart ) THEN ! Restart from a file 87 105 ! ! ------------------- … … 100 118 ! 101 119 ssh(:,:,Kbb) = 0._wp ! set the ocean at rest 120 uu (:,:,:,Kbb) = 0._wp 121 vv (:,:,:,Kbb) = 0._wp 122 ! 102 123 IF( ll_wd ) THEN 103 124 ssh(:,:,Kbb) = -ssh_ref ! Added in 30 here for bathy that adds 30 as Iterative test CEOD … … 111 132 END_2D 112 133 ENDIF 113 uu (:,:,:,Kbb) = 0._wp 114 vv (:,:,:,Kbb) = 0._wp 115 ! 134 ! 116 135 ELSE ! user defined initial T and S 117 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 136 DO jk = 1, jpk 137 zgdept(:,:,jk) = gdept(:,:,jk,Kbb) 138 END DO 139 CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 118 140 ENDIF 119 141 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones … … 121 143 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 122 144 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 123 hdiv(:,:,jpk) = 0._wp ! bottom divergence set one for 0 to zero at jpk level124 CALL div_hor( 0, Kbb, Kmm ) ! compute interior hdiv value125 !!gm hdiv(:,:,:) = 0._wp126 145 127 146 !!gm POTENTIAL BUG : 128 147 !!gm ISSUE : if ssh(:,:,Kbb) /= 0 then, in non linear free surface, the e3._n, e3._b should be recomputed 129 !! as well as gdept and gdepw.... !!!!!148 !! as well as gdept_ and gdepw_.... !!!!! 130 149 !! ===>>>> probably a call to domvvl initialisation here.... 131 150 … … 151 170 ! 152 171 ENDIF 172 #if defined key_agrif 173 ENDIF 174 #endif 153 175 ! 154 176 ! Initialize "now" and "before" barotropic velocities: -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/sbcfwb.F90
r13632 r13877 17 17 USE dom_oce ! ocean space and time domain 18 18 USE sbc_oce ! surface ocean boundary condition 19 USE isf_oce 19 USE isf_oce , ONLY : fwfisf_cav, fwfisf_par, ln_isfcpl, ln_isfcpl_cons, risfcpl_cons_ssh ! ice shelf melting contribution 20 20 USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass 21 21 USE phycst ! physical constants … … 71 71 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmsk_tospread, zerp_cor ! - - 72 72 REAL(wp) ,DIMENSION(1) :: z_fwfprv 73 COMPLEX( wp),DIMENSION(1) :: y_fwfnow73 COMPLEX(dp),DIMENSION(1) :: y_fwfnow 74 74 !!---------------------------------------------------------------------- 75 75 ! … … 207 207 !!gm ===>>>> lbc_lnk should be useless as all the computation is done over the whole domain ! 208 208 #if defined key_mpi3 209 CALL lbc_lnk_nc_ multi( 'sbcfwb', zerp_cor, 'T', 1.)209 CALL lbc_lnk_nc_lnk( 'sbcfwb', zerp_cor, 'T', 1.0_wp ) 210 210 #else 211 CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1. )211 CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1.0_wp ) 212 212 #endif 213 213 ! -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP/EXPREF/namelist_cfg
r13632 r13877 437 437 / 438 438 !----------------------------------------------------------------------- 439 !-----------------------------------------------------------------------440 /441 !-----------------------------------------------------------------------442 439 &namhsb ! Heat and salt budgets (default: OFF) 443 440 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/STATION_ASF/MY_SRC/stpctl.F90
r13632 r13877 65 65 REAL(wp) :: zzz ! local real 66 66 REAL(wp), DIMENSION(4) :: zmax, zmaxlocal 67 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 67 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce 68 68 LOGICAL, DIMENSION(jpi,jpj) :: llmsk 69 69 CHARACTER(len=20) :: clname … … 115 115 llmsk(:,Nje1: jpj) = .FALSE. 116 116 ! 117 llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp ! test only the inner domain 118 IF( COUNT( llmsk(:,:) ) > 0 ) THEN ! avoid huge values sent back for land processors... 119 zmax(1) = MAXVAL( taum(:,:) , mask = llmsk ) ! max wind stress module 120 zmax(2) = MAXVAL( ABS( qns(:,:) ) , mask = llmsk ) ! max non-solar heat flux 121 zmax(3) = MAXVAL( ABS( emp(:,:) ) , mask = llmsk ) ! max E-P 122 ELSE 123 IF( ll_colruns ) THEN ! default value: must not be kept when calling mpp_max -> must be as small as possible 124 zmax(1:3) = -HUGE(1._wp) 125 ELSE ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 126 zmax(1:3) = 0._wp 127 ENDIF 128 ENDIF 129 zmax(4) = REAL( nstop, wp ) ! stop indicator 117 llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp ! test only the inner domain 118 ! 119 ll_0oce = .NOT. ANY( llmsk(:,:) ) ! no ocean point in the inner domain? 120 ! 121 zmax(1) = MAXVAL( taum(:,:) , mask = llmsk ) ! max wind stress module 122 zmax(2) = MAXVAL( ABS( qns(:,:) ), mask = llmsk ) ! max non-solar heat flux 123 zmax(3) = MAXVAL( ABS( emp(:,:) ), mask = llmsk ) ! max E-P 124 zmax(4) = REAL( nstop, wp ) ! stop indicator 125 ! 130 126 ! !== get global extrema ==! 131 127 ! !== done by all processes if writting run.stat ==! … … 134 130 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 135 131 nstop = NINT( zmax(4) ) ! update nstop indicator (now sheared among all local domains) 136 ENDIF 132 ELSE 133 ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow. 134 IF( ll_0oce ) zmax(1:3) = 0._wp ! default "valid" values... 135 ENDIF 136 ! !== error handling ==! 137 137 ! !== write "run.stat" files ==! 138 138 ! !== done only by 1st subdomain at writting timestep ==! 139 139 IF( ll_wrtruns ) THEN 140 140 WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3) 141 istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )142 istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )143 istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/ zmax(3)/), (/kt/), (/1/) )141 DO ji = 1, 3 142 istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) ) 143 END DO 144 144 IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 145 145 END IF
Note: See TracChangeset
for help on using the changeset viewer.