Changeset 13220
- Timestamp:
- 2020-07-02T13:02:36+02:00 (5 years ago)
- Location:
- NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation
- Files:
-
- 56 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation
- Property svn:externals
-
old new 2 2 ^/utils/build/makenemo@HEAD makenemo 3 3 ^/utils/build/mk@HEAD mk 4 ^/utils/tools @HEADtools5 ^/vendors/AGRIF/dev @HEADext/AGRIF4 ^/utils/tools/@HEAD tools 5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg
r13135 r13220 95 95 ! ! bulk algorithm : 96 96 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 97 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003)98 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013)99 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31)100 !101 rn_zqt = 10. ! Air temperature & humidity reference height (m)102 rn_zu = 10. ! Wind vector reference height (m)103 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012)104 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015)105 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow)106 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.)107 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to108 ! ! calculate the wind stress (0.=absolute or 1.=relative winds)109 ln_skin_cs = .false. ! use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB110 ln_skin_wl = .false. ! use the warm-layer " " "111 !112 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true.113 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true.114 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true.115 97 ! 116 98 cn_dir = './' ! root directory for the bulk data location -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
r13135 r13220 92 92 ! ! bulk algorithm : 93 93 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 94 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003)95 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013)96 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31)97 !98 rn_zqt = 10. ! Air temperature & humidity reference height (m)99 rn_zu = 10. ! Wind vector reference height (m)100 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012)101 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015)102 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow)103 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.)104 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to105 ! ! calculate the wind stress (0.=absolute or 1.=relative winds)106 ln_skin_cs = .false. ! use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB107 ln_skin_wl = .false. ! use the warm-layer " " "108 !109 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true.110 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true.111 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true.112 94 ! 113 95 cn_dir = './' ! root directory for the bulk data location … … 176 158 !----------------------------------------------------------------------- 177 159 ln_spc_dyn = .true. ! use 0 as special value for dynamics 178 rn_sponge_tra = 1440. ! coefficient for tracer sponge layer [m2/s]179 rn_sponge_dyn = 1440. ! coefficient for dynamics sponge layer [m2/s]180 160 ln_chk_bathy = .true. ! =T check the parent bathymetry 181 161 / -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
r13135 r13220 158 158 !----------------------------------------------------------------------- 159 159 ln_spc_dyn = .true. ! use 0 as special value for dynamics 160 rn_sponge_tra = 480. ! coefficient for tracer sponge layer [m2/s]161 rn_sponge_dyn = 480. ! coefficient for dynamics sponge layer [m2/s]162 160 ln_chk_bathy = .true. ! =T check the parent bathymetry 163 161 / -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
r13135 r13220 95 95 ! ! bulk algorithm : 96 96 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 97 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003)98 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013)99 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31)100 !101 rn_zqt = 10. ! Air temperature & humidity reference height (m)102 rn_zu = 10. ! Wind vector reference height (m)103 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012)104 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015)105 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow)106 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.)107 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to108 ! ! calculate the wind stress (0.=absolute or 1.=relative winds)109 ln_skin_cs = .false. ! use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB110 ln_skin_wl = .false. ! use the warm-layer " " "111 !112 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true.113 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true.114 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true.115 97 ! 116 98 cn_dir = './' ! root directory for the bulk data location -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml
r12063 r13220 56 56 <field field_ref="t_abl" /> 57 57 <field field_ref="q_abl" /> 58 <field field_ref="uvz1_abl" /> 59 <field field_ref="tz1_abl" /> 60 <field field_ref="qz1_abl" /> 61 <field field_ref="uvz1_dta" /> 62 <field field_ref="tz1_dta" /> 63 <field field_ref="qz1_dta" /> 58 64 <field field_ref="pblh" /> 59 65 </file> -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg
r13135 r13220 110 110 ! ! bulk algorithm : 111 111 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 112 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003)113 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013)114 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31)115 rn_zqt = 10. ! Air temperature & humidity reference height (m)116 rn_zu = 10. ! Wind vector reference height (m)117 !118 ! Skin is ONLY available in ECMWF and COARE algorithms:119 ln_skin_cs = .false. ! use the cool-skin parameterization => set nn_fsbc=1 and ln_dm2dc=.true.!120 ln_skin_wl = .false. ! use the warm-layer " => set nn_fsbc=1 and ln_dm2dc=.true.!121 !122 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true.123 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true.124 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true.125 112 ! 126 113 cn_dir = './' ! root directory for the bulk data location … … 132 119 sn_tair = 'tair_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d', 24., 'tair' , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bilinear' , '' , '' 133 120 sn_humi = 'humi_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d', 24., 'humi' , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bilinear' , '' , '' 134 sn_hpgi = 'uhpg_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d', 24., 'uhpg' , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic' , 'UG' , ''135 sn_hpgj = 'vhpg_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d', 24., 'vhpg' , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic' , 'VG' , ''136 137 121 sn_qsr = 'ncar_rad.15JUNE2009_fill' , 24., 'SWDN_MOD', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 138 122 sn_qlw = 'ncar_rad.15JUNE2009_fill' , 24., 'LWDN_MOD', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' … … 140 124 sn_snow = 'ncar_precip.15JUNE2009_fill' , -1., 'SNOW' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 141 125 sn_slp = 'slp.15JUNE2009_fill' , 6., 'SLP' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 126 sn_hpgi = 'uhpg_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d', 24., 'uhpg' , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic' , 'UG' , '' 127 sn_hpgj = 'vhpg_drwnlnd_ERAI_L25Z10_GLOBAL_F128R_ana1d', 24., 'vhpg' , .false. , .false. , 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic' , 'VG' , '' 142 128 / 143 129 -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
r13135 r13220 121 121 / 122 122 !----------------------------------------------------------------------- 123 &namsbc_abl ! Atmospheric Boundary Layer formulation (ln_abl = T) 124 !----------------------------------------------------------------------- 125 / 126 !----------------------------------------------------------------------- 123 127 &namtra_qsr ! penetrative solar radiation (ln_traqsr =T) 124 128 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg
r13135 r13220 68 68 ! ! bulk algorithm : 69 69 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 70 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003)71 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013)72 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31)73 !74 rn_zqt = 10. ! Air temperature & humidity reference height (m)75 rn_zu = 10. ! Wind vector reference height (m)76 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012)77 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015)78 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow)79 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.)80 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to81 ! ! calculate the wind stress (0.=absolute or 1.=relative winds)82 ln_skin_cs = .false. ! use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB83 ln_skin_wl = .false. ! use the warm-layer " " "84 !85 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true.86 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true.87 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true.88 70 ! 89 71 cn_dir = './' ! root directory for the bulk data location -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/SHARED/field_def_nemo-oce.xml
r13135 r13220 455 455 <field id="t_dta" long_name="DTA potential temperature" standard_name="dta_theta" unit="K" /> 456 456 <field id="q_dta" long_name="DTA specific humidity" standard_name="dta_qspe" unit="kg/kg" /> 457 <field id="coeft" long_name="ABL nudging coefficient" standard_name="coeft" unit="" /> 457 <field id="u_geo" long_name="GEO i-horizontal velocity" standard_name="geo_x_velocity" unit="m/s" /> 458 <field id="v_geo" long_name="GEO j-horizontal velocity" standard_name="geo_y_velocity" unit="m/s" /> 458 459 <field id="tke_abl" long_name="ABL turbulent kinetic energy" standard_name="abl_tke" unit="m2/s2" /> 459 460 <field id="avm_abl" long_name="ABL turbulent viscosity" standard_name="abl_avm" unit="m2/s" /> 460 461 <field id="avt_abl" long_name="ABL turbulent diffusivity" standard_name="abl_avt" unit="m2/s" /> 461 <field id="mxl_abl" long_name="ABL mixing length" standard_name="abl_mxl" unit="m" /> 462 <field id="mxlm_abl" long_name="ABL master mixing length" standard_name="abl_mxlm" unit="m" /> 463 <field id="mxld_abl" long_name="ABL dissipative mixing length" standard_name="abl_mxld" unit="m" /> 462 464 </field_group> 463 465 -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/SHARED/namelist_ref
r13135 r13220 268 268 ln_Cd_L12 = .false. ! air-ice drags = F(ice conc.) (Lupkes et al. 2012) 269 269 ln_Cd_L15 = .false. ! air-ice drags = F(ice conc.) (Lupkes et al. 2015) 270 ! ! - module of the mean stress" data 270 ln_crt_fbk = .false. ! Add surface current feedback to the wind stress (Renault et al. 2020, doi: 10.1029/2019MS001715) 271 rn_stau_a = -2.9e-3 ! Alpha from eq. 10: Stau = Alpha * Wnd + Beta 272 rn_stau_b = 8.0e-3 ! Beta 271 273 rn_pfac = 1. ! multipl. factor for precipitation (total & snow) 272 274 rn_efac = 1. ! multipl. factor for evaporation (0. or 1.) 273 rn_vfac = 0. ! multipl. factor for ocean & ice velocity274 ! ! used to calculate the wind stress275 ! ! (0. => absolute or 1. => relative winds)276 275 ln_skin_cs = .false. ! use the cool-skin parameterization 277 276 ln_skin_wl = .false. ! use the warm-layer parameterization … … 280 279 ln_humi_dpt = .false. ! humidity "sn_humi" is dew-point temperature [K] 281 280 ln_humi_rlh = .false. ! humidity "sn_humi" is relative humidity [%] 281 ln_tpot = .true. !!GS: compute potential temperature or not 282 282 ! 283 283 cn_dir = './' ! root directory for the bulk data location … … 291 291 sn_tair = 't_10.15JUNE2009_fill' , 6. , 'T_10_MOD', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 292 292 sn_humi = 'q_10.15JUNE2009_fill' , 6. , 'Q_10_MOD', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 293 sn_hpgi = 'NOT USED' , 24. , 'uhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG' , ''294 sn_hpgj = 'NOT USED' , 24. , 'vhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG' , ''295 293 sn_prec = 'ncar_precip.15JUNE2009_fill', -1. , 'PRC_MOD1', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 296 294 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 297 295 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 296 sn_uoatm = 'NOT USED' , 6. , 'UOATM' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , 'Uoceatm', '' 297 sn_voatm = 'NOT USED' , 6. , 'VOATM' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , 'Voceatm', '' 298 sn_hpgi = 'NOT USED' , 24. , 'uhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG' , '' 299 sn_hpgj = 'NOT USED' , 24. , 'vhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG' , '' 298 300 / 299 301 !----------------------------------------------------------------------- … … 308 310 cn_ablrst_outdir = "." ! directory to write output abl restarts 309 311 312 ln_rstart_abl = .false. 310 313 ln_hpgls_frc = .false. 311 314 ln_geos_winds = .false. 312 nn_dyn_restore = 2 ! restoring option for dynamical ABL variables: = 0 no restoring 315 ln_smth_pblh = .false. 316 nn_dyn_restore = 0 ! restoring option for dynamical ABL variables: = 0 no restoring 313 317 ! = 1 equatorial restoring 314 318 ! = 2 global restoring 315 rn_ldyn_min = 4.5 ! magnitude of the nudging on ABL dynamics at the bottom of the ABL [hour]316 rn_ldyn_max = 1.5 ! magnitude of the nudging on ABL dynamics at the top of the ABL [hour]317 rn_ltra_min = 4.5 ! magnitude of the nudging on ABL tracers at the bottom of the ABL [hour]318 rn_ltra_max = 1.5 ! magnitude of the nudging on ABL tracers at the top of the ABL [hour]319 rn_ldyn_min = 4.5 ! dynamics nudging magnitude inside the ABL [hour] (~3 rn_Dt) 320 rn_ldyn_max = 1.5 ! dynamics nudging magnitude above the ABL [hour] (~1 rn_Dt) 321 rn_ltra_min = 4.5 ! tracers nudging magnitude inside the ABL [hour] (~3 rn_Dt) 322 rn_ltra_max = 1.5 ! tracers nudging magnitude above the ABL [hour] (~1 rn_Dt) 319 323 nn_amxl = 0 ! mixing length: = 0 Deardorff 80 length-scale 320 324 ! = 1 length-scale based on the distance to the PBL height 321 325 ! = 2 Bougeault & Lacarrere 89 length-scale 322 rn_Cm = 0.0667 ! 0.126 in MesoNH 323 rn_Ct = 0.1667 ! 0.143 in MesoNH 324 rn_Ce = 0.4 ! 0.4 in MesoNH 325 rn_Ceps = 0.7 ! 0.85 in MesoNH 326 rn_Rod = 0.15 ! c0 in RMCA17 mixing length formulation (not yet implemented) 327 rn_Ric = 0.139 ! Critical Richardson number (to compute PBL height and diffusivities) 326 ! CBR00 ! CCH02 ! MesoNH ! 327 rn_Cm = 0.0667 ! 0.0667 ! 0.1260 ! 0.1260 ! 328 rn_Ct = 0.1667 ! 0.1667 ! 0.1430 ! 0.1430 ! 329 rn_Ce = 0.40 ! 0.40 ! 0.34 ! 0.40 ! 330 rn_Ceps = 0.700 ! 0.700 ! 0.845 ! 0.850 ! 331 rn_Ric = 0.139 ! 0.139 ! 0.143 ! ? ! Critical Richardson number (to compute PBL height and diffusivities) 332 rn_Rod = 0.15 ! c0 in RMCA17 mixing length formulation (not yet implemented) 328 333 / 329 334 !----------------------------------------------------------------------- … … 638 643 &namagrif ! AGRIF zoom ("key_agrif") 639 644 !----------------------------------------------------------------------- 640 ln_agrif_2way = .true. ! activate two way nesting 641 ln_spc_dyn = .true. ! use 0 as special value for dynamics 642 rn_sponge_tra = 2880. ! coefficient for tracer sponge layer [m2/s] 643 rn_sponge_dyn = 2880. ! coefficient for dynamics sponge layer [m2/s] 644 rn_trelax_tra = 0.01 ! inverse of relaxation time (in steps) for tracers [] 645 rn_trelax_dyn = 0.01 ! inverse of relaxation time (in steps) for dynamics [] 646 ln_chk_bathy = .false. ! =T check the parent bathymetry 645 ln_agrif_2way = .true. ! activate two way nesting 646 ln_init_chfrpar = .false. ! initialize child grids from parent 647 ln_spc_dyn = .true. ! use 0 as special value for dynamics 648 rn_sponge_tra = 0.002 ! coefficient for tracer sponge layer [] 649 rn_sponge_dyn = 0.002 ! coefficient for dynamics sponge layer [] 650 rn_trelax_tra = 0.01 ! inverse of relaxation time (in steps) for tracers [] 651 rn_trelax_dyn = 0.01 ! inverse of relaxation time (in steps) for dynamics [] 652 ln_chk_bathy = .false. ! =T check the parent bathymetry 647 653 / 648 654 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/cfgs/WED025/EXPREF/namelist_cfg
r13135 r13220 138 138 !----------------------------------------------------------------------- 139 139 ! ! bulk algorithm : 140 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008)140 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 141 141 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 142 142 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013) 143 143 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 45r1) 144 144 ! 145 145 cn_dir = './' ! root directory for the bulk data location 146 146 !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/doc/latex/NEMO/subfiles/chap_SBC.tex
r12377 r13220 832 832 Solid precipitation & snow & $Kg.m^{-2}.s^{-1}$ & T \\ 833 833 \hline 834 Mean sea-level pressure & slp & $ hPa$ & T \\834 Mean sea-level pressure & slp & $Pa$ & T \\ 835 835 \hline 836 836 \end{tabular} -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ABL/abl.F90
r12489 r13220 29 29 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: avm_abl !: turbulent viscosity [m2/s] 30 30 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: avt_abl !: turbulent diffusivity [m2/s] 31 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: mxl_abl !: mixing length [m] 31 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: mxld_abl !: dissipative mixing length [m] 32 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: mxlm_abl !: master mixing length [m] 32 33 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: tke_abl !: turbulent kinetic energy [m2/s2] 33 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: fft_abl !: Coriolis parameter [1/s] … … 55 56 !!---------------------------------------------------------------------- 56 57 ! 57 ALLOCATE( u_abl (1:jpi,1:jpj,1:jpka,jptime), & 58 & v_abl (1:jpi,1:jpj,1:jpka,jptime), & 59 & tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), & 60 & avm_abl(1:jpi,1:jpj,1:jpka), & 61 & avt_abl(1:jpi,1:jpj,1:jpka), & 62 & mxl_abl(1:jpi,1:jpj,1:jpka), & 63 & tke_abl(1:jpi,1:jpj,1:jpka,jptime), & 64 & fft_abl(1:jpi,1:jpj), & 65 & pblh (1:jpi,1:jpj), & 66 & msk_abl(1:jpi,1:jpj), & 67 & rest_eq(1:jpi,1:jpj), & 68 & e3t_abl(1:jpka), e3w_abl(1:jpka), ght_abl(1:jpka), ghw_abl(1:jpka), STAT=ierr ) 58 ALLOCATE( u_abl (1:jpi,1:jpj,1:jpka,jptime ), & 59 & v_abl (1:jpi,1:jpj,1:jpka,jptime ), & 60 & tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), & 61 & tke_abl (1:jpi,1:jpj,1:jpka,jptime ), & 62 & avm_abl (1:jpi,1:jpj,1:jpka ), & 63 & avt_abl (1:jpi,1:jpj,1:jpka ), & 64 & mxld_abl(1:jpi,1:jpj,1:jpka ), & 65 & mxlm_abl(1:jpi,1:jpj,1:jpka ), & 66 & fft_abl (1:jpi,1:jpj ), & 67 & pblh (1:jpi,1:jpj ), & 68 & msk_abl (1:jpi,1:jpj ), & 69 & rest_eq (1:jpi,1:jpj ), & 70 & e3t_abl (1:jpka), e3w_abl(1:jpka) , & 71 & ght_abl (1:jpka), ghw_abl(1:jpka) , STAT=ierr ) 69 72 ! 70 73 abl_alloc = ierr -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ABL/ablmod.F90
r13135 r13220 2 2 !!====================================================================== 3 3 !! *** MODULE ablmod *** 4 !! Surface module : ABL computation to provide atmospheric data 4 !! Surface module : ABL computation to provide atmospheric data 5 5 !! for surface fluxes computation 6 6 !!====================================================================== 7 7 !! History : 3.6 ! 2019-03 (F. Lemarié & G. Samson) Original code 8 8 !!---------------------------------------------------------------------- 9 9 10 10 !!---------------------------------------------------------------------- 11 11 !! abl_stp : ABL single column model … … 16 16 17 17 USE phycst ! physical constants 18 USE dom_oce, ONLY : tmask 18 USE dom_oce, ONLY : tmask 19 19 USE sbc_oce, ONLY : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa 20 USE sbcblk ! use rn_ ?fac20 USE sbcblk ! use rn_efac, cdn_oce 21 21 USE sbcblk_phy ! use some physical constants for flux computation 22 22 ! … … 30 30 31 31 PUBLIC abl_stp ! called by sbcabl.F90 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2, zrough 33 33 !! * Substitutions 34 34 # include "do_loop_substitute.h90" … … 38 38 39 39 !=================================================================================================== 40 SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &! in41 & pu_dta, pv_dta, pt_dta, pq_dta, & 40 SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, & ! in 41 & pu_dta, pv_dta, pt_dta, pq_dta, & 42 42 & pslp_dta, pgu_dta, pgv_dta, & 43 & pcd_du, psen, pevp, & ! in/out 44 & pwndm, ptaui, ptauj, ptaum & 45 #if defined key_si3 46 & , ptm_su,pssu_ice,pssv_ice,pssq_ice,pcd_du_ice & 47 & , psen_ice, pevp_ice, pwndm_ice, pfrac_oce & 48 & , ptaui_ice, ptauj_ice & 49 #endif 50 & ) 43 & pcd_du, psen, pevp, & ! in/out 44 & pwndm, ptaui, ptauj, ptaum & 45 #if defined key_si3 46 & , ptm_su, pssu_ice, pssv_ice & 47 & , pssq_ice, pcd_du_ice, psen_ice & 48 & , pevp_ice, pwndm_ice, pfrac_oce & 49 & , ptaui_ice, ptauj_ice & 50 #endif 51 & ) 51 52 !--------------------------------------------------------------------------------------------------- 52 53 … … 54 55 !! *** ROUTINE abl_stp *** 55 56 !! 56 !! ** Purpose : Time-integration of the ABL model 57 !! ** Purpose : Time-integration of the ABL model 57 58 !! 58 !! ** Method : Compute atmospheric variables : vertical turbulence 59 !! ** Method : Compute atmospheric variables : vertical turbulence 59 60 !! + Coriolis term + newtonian relaxation 60 !! 61 !! 61 62 !! ** Action : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh 62 63 !! - Advance tracers to time n+1 (Euler backward scheme) … … 70 71 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psst ! sea-surface temperature [Celsius] 71 72 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu ! sea-surface u (U-point) 72 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv ! sea-surface v (V-point) 73 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv ! sea-surface v (V-point) 73 74 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq ! sea-surface humidity 74 75 REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pu_dta ! large-scale windi … … 82 83 REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: psen ! Ch x Du 83 84 REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pevp ! Ce x Du 84 REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! ||uwnd|| 85 REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! ||uwnd|| 85 86 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui ! taux 86 87 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj ! tauy 87 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaum ! ||tau|| 88 ! 89 #if defined key_si3 88 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaum ! ||tau|| 89 ! 90 #if defined key_si3 90 91 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptm_su ! ice-surface temperature [K] 91 92 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu_ice ! ice-surface u (U-point) 92 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv_ice ! ice-surface v (V-point) 93 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq_ice ! ice-surface humidity 93 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv_ice ! ice-surface v (V-point) 94 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq_ice ! ice-surface humidity 94 95 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du_ice ! Cd x Du over ice (T-point) 95 96 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psen_ice ! Ch x Du over ice (T-point) 96 97 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pevp_ice ! Ce x Du over ice (T-point) 97 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndm_ice ! ||uwnd - uice|| 98 !REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pfrac_oce !!GS: out useless ? 99 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pfrac_oce ! 98 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndm_ice ! ||uwnd - uice|| 99 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pfrac_oce ! ocean fraction 100 100 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui_ice ! ice-surface taux stress (U-point) 101 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point) 102 #endif 103 ! 104 REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j 105 REAL(wp), DIMENSION(1:jpi,2:jpka ) :: zCF 106 REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: z_cft !--FL--to be removed after the test phase 107 ! 108 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a 109 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b 110 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c 101 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point) 102 #endif 103 ! 104 REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j 105 REAL(wp), DIMENSION(1:jpi ,2:jpka) :: zCF 106 ! 107 REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_a 108 REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_b 109 REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_c 111 110 ! 112 111 INTEGER :: ji, jj, jk, jtra, jbak ! dummy loop indices 113 112 REAL(wp) :: zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2 114 113 REAL(wp) :: zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce 115 ! 116 !!--------------------------------------------------------------------- 114 LOGICAL :: SemiImp_Cor = .TRUE. 115 ! 116 !!--------------------------------------------------------------------- 117 117 ! 118 118 IF(lwp .AND. kt == nit000) THEN ! control print … … 120 120 WRITE(numout,*) 'abl_stp : ABL time stepping' 121 121 WRITE(numout,*) '~~~~~~' 122 ENDIF 122 ENDIF 123 123 ! 124 124 IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) ) 125 !! Compute ustar squared as Cd || Uatm-Uoce ||^2 126 !! needed for surface boundary condition of TKE 125 IF( kt == nit000 ) ALLOCATE ( zrough( 1:jpi, 1:jpj ) ) 126 !! Compute ustar squared as Cd || Uatm-Uoce ||^2 127 !! needed for surface boundary condition of TKE 127 128 !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk) 128 129 DO_2D_11_11 129 130 zzoce = pCd_du (ji,jj) * pwndm (ji,jj) 130 131 #if defined key_si3 131 zzice = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 132 ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 132 zzice = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 133 ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 133 134 #else 134 ustar2(ji,jj) = zzoce 135 ustar2(ji,jj) = zzoce 135 136 #endif 137 zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(ji,jj), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce 136 138 END_2D 137 139 ! … … 140 142 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 141 143 142 CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj) 143 144 CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj) 145 144 146 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 145 147 ! ! 2 *** Advance tracers to time n+1 146 148 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 147 149 148 150 !------------- 149 151 DO jj = 1, jpj ! outer loop !--> tq_abl computed on (1:jpi) x (1:jpj) 150 !------------- 151 ! Compute matrix elements for interior points 152 !------------- 153 ! Compute matrix elements for interior points 152 154 DO jk = 3, jpkam1 153 155 DO ji = 1, jpi ! vector opt. 154 z_elem_a( ji, jk) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal155 z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal156 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal157 END DO 158 END DO 159 ! Boundary conditions 160 DO ji = 1, jpi ! vector opt. 161 ! Neumann at the bottom 162 z_elem_a( ji, 2) = 0._wp163 z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 )156 z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 157 z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 158 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 159 END DO 160 END DO 161 ! Boundary conditions 162 DO ji = 1, jpi ! vector opt. 163 ! Neumann at the bottom 164 z_elem_a( ji, 2 ) = 0._wp 165 z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 ) 164 166 ! Homogeneous Neumann at the top 165 z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )166 z_elem_c( ji, jpka) = 0._wp167 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji,jpka )168 END DO 167 z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 168 z_elem_c( ji, jpka ) = 0._wp 169 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 170 END DO 169 171 170 172 DO jtra = 1,jptq ! loop on active tracers 171 173 172 174 DO jk = 3, jpkam1 173 DO ji = 1,jpi 174 tq_abl ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl ( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side 175 !DO ji = 2, jpim1 176 DO ji = 1,jpi !!GS: to be checked if needed 177 tq_abl( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side 175 178 END DO 176 179 END DO 177 180 178 181 IF(jtra == jp_ta) THEN 179 DO ji = 1,jpi ! boundary conditions for temperature180 zztmp1 = psen(ji, jj) 181 zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 ) 182 #if defined key_si3 183 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 182 DO ji = 1,jpi ! surface boundary condition for temperature 183 zztmp1 = psen(ji, jj) 184 zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 ) 185 #if defined key_si3 186 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 184 187 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj) 185 #endif 186 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1187 tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2 , nt_n, jtra ) + rDt_abl * zztmp2188 #endif 189 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 190 tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2 188 191 tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra ) 189 END DO 190 ELSE 191 DO ji = 1,jpi ! boundary conditions for humidity192 zztmp1 = pevp(ji, jj) 193 zztmp2 = pevp(ji, jj) * pssq(ji, jj) 194 #if defined key_si3 195 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 196 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj) 197 #endif 198 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2) + rDt_abl * zztmp1199 tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2 , nt_n, jtra ) + rDt_abl * zztmp2192 END DO 193 ELSE ! jp_qa 194 DO ji = 1,jpi ! surface boundary condition for humidity 195 zztmp1 = pevp(ji, jj) 196 zztmp2 = pevp(ji, jj) * pssq(ji, jj) 197 #if defined key_si3 198 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 199 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj) 200 #endif 201 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 202 tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2 200 203 tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra ) 201 END DO 204 END DO 202 205 END IF 203 206 !! 204 207 !! Matrix inversion 205 208 !! ---------------------------------------------------------- 206 DO ji = 1,jpi 207 zcff = 1._wp / z_elem_b( ji, 2 )208 zCF ( ji, 2) = - zcff * z_elem_c( ji, 2 )209 tq_abl( ji,jj,2,nt_a,jtra) = zcff * tq_abl(ji,jj,2,nt_a,jtra)210 END DO 211 212 DO jk = 3, jpka 213 DO ji = 1,jpi 214 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )209 DO ji = 1,jpi 210 zcff = 1._wp / z_elem_b( ji, 2 ) 211 zCF ( ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 212 tq_abl( ji, jj, 2, nt_a, jtra ) = zcff * tq_abl( ji, jj, 2, nt_a, jtra ) 213 END DO 214 215 DO jk = 3, jpka 216 DO ji = 1,jpi 217 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF( ji, jk-1 ) ) 215 218 zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) 216 219 tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk ,nt_a,jtra) & 217 & - z_elem_a(ji, jk) *tq_abl(ji,jj,jk-1,nt_a,jtra) )218 END DO 219 END DO 220 !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ... 221 DO jk = jpkam1,2,-1 222 DO ji = 1,jpi 220 & - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) ) 221 END DO 222 END DO 223 !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ... 224 DO jk = jpkam1,2,-1 225 DO ji = 1,jpi 223 226 tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) + & 224 227 & zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra) 225 228 END DO 226 229 END DO 227 228 END DO !<-- loop on tracers 229 !! 230 !------------- 231 END DO ! end outer loop 232 !------------- 233 234 230 231 END DO !<-- loop on tracers 232 !! 233 !------------- 234 END DO ! end outer loop 235 !------------- 236 235 237 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 236 238 ! ! 3 *** Compute Coriolis term with geostrophic guide 237 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 238 !------------- 239 DO jk = 2, jpka ! outer loop 240 !------------- 241 ! 242 ! Advance u_abl & v_abl to time n+1 243 DO_2D_11_11 244 zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl ) ! (f dt)**2 239 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 240 IF( SemiImp_Cor ) THEN 241 242 !------------- 243 DO jk = 2, jpka ! outer loop 244 !------------- 245 ! 246 ! Advance u_abl & v_abl to time n+1 247 DO_2D_11_11 248 zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl ) ! (f dt)**2 249 250 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 251 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * u_abl( ji, jj, jk, nt_n ) & 252 & + rDt_abl * fft_abl(ji, jj) * v_abl( ji, jj, jk, nt_n ) ) & 253 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 245 254 246 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 247 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n ) & 248 & + rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj , jk, nt_n ) ) & 249 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 250 251 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 252 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n ) & 253 & - rDt_abl * fft_abl(ji, jj) * u_abl ( ji , jj, jk, nt_n ) ) & 254 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 255 END_2D 256 ! 257 !------------- 258 END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 259 !------------- 260 ! 261 IF( ln_geos_winds ) THEN 262 DO jj = 1, jpj ! outer loop 263 DO jk = 1, jpka 264 DO ji = 1, jpi 265 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) & 266 & - rDt_abl * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk) 267 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) & 268 & + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj ) * pgu_dta(ji ,jj ,jk) 269 END DO 270 END DO 271 END DO 272 END IF 273 !------------- 274 ! 275 IF( ln_hpgls_frc ) THEN 276 DO jj = 1, jpj ! outer loop 277 DO jk = 1, jpka 278 DO ji = 1, jpi 279 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 280 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 255 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 256 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * v_abl( ji, jj, jk, nt_n ) & 257 & - rDt_abl * fft_abl(ji, jj) * u_abl( ji, jj, jk, nt_n ) ) & 258 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 259 END_2D 260 ! 261 !------------- 262 END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 263 !------------- 264 ! 265 IF( ln_geos_winds ) THEN 266 DO jj = 1, jpj ! outer loop 267 DO jk = 1, jpka 268 DO ji = 1, jpi 269 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) & 270 & - rDt_abl * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk) 271 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) & 272 & + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj ) * pgu_dta(ji ,jj ,jk) 273 END DO 274 END DO 275 END DO 276 END IF 277 ! 278 IF( ln_hpgls_frc ) THEN 279 DO jj = 1, jpj ! outer loop 280 DO jk = 1, jpka 281 DO ji = 1, jpi 282 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 283 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 284 ENDDO 281 285 ENDDO 282 286 ENDDO 283 ENDDO 284 END IF 285 286 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 287 ! ! 4 *** Advance u,v to time n+1 288 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 289 ! 290 ! Vertical diffusion for u_abl 287 END IF 288 289 ELSE ! SemiImp_Cor = .FALSE. 290 291 IF( ln_geos_winds ) THEN 292 293 !------------- 294 DO jk = 2, jpka ! outer loop 295 !------------- 296 ! 297 IF( MOD( kt, 2 ) == 0 ) then 298 ! Advance u_abl & v_abl to time n+1 299 DO jj = 1, jpj 300 DO ji = 1, jpi 301 zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_n ) - pgv_dta(ji ,jj ,jk) ) 302 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff 303 zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_a ) - pgu_dta(ji ,jj ,jk) ) 304 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff ) 305 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * u_abl( ji, jj, jk, nt_a ) 306 END DO 307 END DO 308 ELSE 309 DO jj = 1, jpj 310 DO ji = 1, jpi 311 zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_n ) - pgu_dta(ji ,jj ,jk) ) 312 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff 313 zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_a ) - pgv_dta(ji ,jj ,jk) ) 314 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff ) 315 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * v_abl( ji, jj, jk, nt_a ) 316 END DO 317 END DO 318 END IF 319 ! 320 !------------- 321 END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 322 !------------- 323 324 ENDIF ! ln_geos_winds 325 326 ENDIF ! SemiImp_Cor 327 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 328 ! ! 4 *** Advance u,v to time n+1 329 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 330 ! 331 ! Vertical diffusion for u_abl 291 332 !------------- 292 333 DO jj = 1, jpj ! outer loop 293 !------------- 334 !------------- 294 335 295 336 DO jk = 3, jpkam1 296 DO ji = 1, jpi 297 z_elem_a( ji, 298 z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal299 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )! diagonal300 END DO 301 END DO 302 303 DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi) 337 DO ji = 1, jpi 338 z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 339 z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 340 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 341 END DO 342 END DO 343 344 DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi) 304 345 !++ Surface boundary condition 305 z_elem_a( ji, 2) = 0._wp306 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )307 ! 308 309 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) 310 #if defined key_si3 346 z_elem_a( ji, 2 ) = 0._wp 347 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 348 ! 349 zztmp1 = pcd_du(ji, jj) 350 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) 351 #if defined key_si3 311 352 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) 312 zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) )313 314 #endif 315 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 353 zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj) ) 354 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 355 #endif 356 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 316 357 u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2 317 318 !++ Top Neumann B.C. 319 !z_elem_a( ji, jpka ) = - 0.5_wp * rDt_abl * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka ) 320 !z_elem_c( ji, jpka ) = 0._wp 321 !z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 322 !++ Top Dirichlet B.C. 323 z_elem_a( ji, jpka ) = 0._wp 324 z_elem_c( ji, jpka ) = 0._wp 325 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 326 u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 327 END DO 358 359 ! idealized test cases only 360 !IF( ln_topbc_neumann ) THEN 361 ! !++ Top Neumann B.C. 362 ! z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 363 ! z_elem_c( ji, jpka ) = 0._wp 364 ! z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 365 ! !u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl ( ji, jj, jpka, nt_a ) 366 !ELSE 367 !++ Top Dirichlet B.C. 368 z_elem_a( ji, jpka ) = 0._wp 369 z_elem_c( ji, jpka ) = 0._wp 370 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 371 u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 372 !ENDIF 373 374 END DO 328 375 !! 329 376 !! Matrix inversion 330 377 !! ---------------------------------------------------------- 331 DO ji = 2, jpi 378 !DO ji = 2, jpi 379 DO ji = 1, jpi !!GS: TBI 332 380 zcff = 1._wp / z_elem_b( ji, 2 ) 333 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 381 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 334 382 u_abl (ji,jj,2,nt_a) = zcff * u_abl(ji,jj,2,nt_a) 335 END DO 336 337 DO jk = 3, jpka 338 DO ji = 2, jpi 383 END DO 384 385 DO jk = 3, jpka 386 DO ji = 2, jpi 339 387 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) ) 340 388 zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) … … 343 391 END DO 344 392 END DO 345 346 DO jk = jpkam1,2,-1 393 394 DO jk = jpkam1,2,-1 347 395 DO ji = 2, jpi 348 396 u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a) 349 397 END DO 350 398 END DO 351 352 !------------- 353 END DO ! end outer loop 354 !------------- 355 356 ! 357 ! Vertical diffusion for v_abl 399 400 !------------- 401 END DO ! end outer loop 402 !------------- 403 404 ! 405 ! Vertical diffusion for v_abl 358 406 !------------- 359 407 DO jj = 2, jpj ! outer loop 360 !------------- 408 !------------- 361 409 ! 362 410 DO jk = 3, jpkam1 363 DO ji = 1, jpi 364 z_elem_a( ji, 365 z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal366 z_elem_b( ji, 367 END DO 368 END DO 369 370 DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 411 DO ji = 1, jpi 412 z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 413 z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 414 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 415 END DO 416 END DO 417 418 DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 371 419 !++ Surface boundary condition 372 z_elem_a( ji, 2) = 0._wp373 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )374 ! 375 376 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) 377 #if defined key_si3 420 z_elem_a( ji, 2 ) = 0._wp 421 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 422 ! 423 zztmp1 = pcd_du(ji, jj) 424 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) 425 #if defined key_si3 378 426 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) 379 zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) )380 381 #endif 382 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 427 zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) ) 428 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 429 #endif 430 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 383 431 v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2 384 !++ Top Neumann B.C. 385 !z_elem_a( ji, jpka ) = -rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 386 !z_elem_c( ji, jpka ) = 0._wp 387 !z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 388 !++ Top Dirichlet B.C. 389 z_elem_a( ji, jpka ) = 0._wp 390 z_elem_c( ji, jpka ) = 0._wp 391 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 392 v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 393 END DO 432 433 ! idealized test cases only 434 !IF( ln_topbc_neumann ) THEN 435 ! !++ Top Neumann B.C. 436 ! z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 437 ! z_elem_c( ji, jpka ) = 0._wp 438 ! z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 439 ! !v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl ( ji, jj, jpka, nt_a ) 440 !ELSE 441 !++ Top Dirichlet B.C. 442 z_elem_a( ji, jpka ) = 0._wp 443 z_elem_c( ji, jpka ) = 0._wp 444 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 445 v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 446 !ENDIF 447 448 END DO 394 449 !! 395 450 !! Matrix inversion 396 451 !! ---------------------------------------------------------- 397 DO ji = 1, jpi 452 DO ji = 1, jpi 398 453 zcff = 1._wp / z_elem_b( ji, 2 ) 399 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 400 v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a ) 401 END DO 402 403 DO jk = 3, jpka 404 DO ji = 1, jpi 454 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 455 v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a ) 456 END DO 457 458 DO jk = 3, jpka 459 DO ji = 1, jpi 405 460 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) ) 406 461 zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) … … 409 464 END DO 410 465 END DO 411 412 DO jk = jpkam1,2,-1 413 DO ji = 1, jpi 466 467 DO jk = jpkam1,2,-1 468 DO ji = 1, jpi 414 469 v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a) 415 470 END DO 416 471 END DO 417 ! 418 !------------- 419 END DO ! end outer loop 420 !------------- 421 422 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 423 ! ! 5 *** Apply nudging on the dynamics and the tracers 424 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 425 z_cft(:,:,:) = 0._wp 426 472 ! 473 !------------- 474 END DO ! end outer loop 475 !------------- 476 477 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 478 ! ! 5 *** Apply nudging on the dynamics and the tracers 479 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 480 427 481 IF( nn_dyn_restore > 0 ) THEN 428 !------------- 482 !------------- 429 483 DO jk = 2, jpka ! outer loop 430 !------------- 484 !------------- 431 485 DO_2D_01_01 432 486 zcff1 = pblh( ji, jj ) 433 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 434 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 487 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 488 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 435 489 zmsk = msk_abl(ji,jj) 436 490 zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2 & 437 491 & + jp_alp1_dyn * zsig + jp_alp0_dyn 438 492 zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points 439 ! rn_Dt = rDt_abl / nn_fsbc 493 ! rn_Dt = rDt_abl / nn_fsbc 440 494 zcff = zcff * rest_eq(ji,jj) 441 z_cft( ji, jj, jk ) = zcff442 495 u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) & 443 & + zcff * pu_dta( ji, jj, jk ) 496 & + zcff * pu_dta( ji, jj, jk ) 444 497 v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * v_abl( ji, jj, jk, nt_a ) & 445 498 & + zcff * pv_dta( ji, jj, jk ) … … 447 500 !------------- 448 501 END DO ! end outer loop 449 !------------- 502 !------------- 450 503 END IF 451 504 452 !------------- 505 !------------- 453 506 DO jk = 2, jpka ! outer loop 454 !------------- 507 !------------- 455 508 DO_2D_11_11 456 509 zcff1 = pblh( ji, jj ) 457 510 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 458 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 511 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 459 512 zmsk = msk_abl(ji,jj) 460 513 zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2 & 461 514 & + jp_alp1_tra * zsig + jp_alp0_tra 462 515 zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points 463 ! rn_Dt = rDt_abl / nn_fsbc 464 !z_cft( ji, jj, jk ) = zcff 516 ! rn_Dt = rDt_abl / nn_fsbc 465 517 tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) & 466 518 & + zcff * pt_dta( ji, jj, jk ) 467 519 468 520 tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa ) & 469 521 & + zcff * pq_dta( ji, jj, jk ) 470 522 471 523 END_2D 472 524 !------------- 473 525 END DO ! end outer loop 474 !------------- 475 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 476 ! ! 6 *** MPI exchanges 477 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 478 ! 479 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1.0_wp, v_abl(:,:,:,nt_a ), 'T', -1.0_wp ) 480 CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1.0_wp, tq_abl(:,:,:,nt_a,jp_qa), 'T', 1.0_wp, kfillmode = jpfillnothing ) ! ++++ this should not be needed... 481 ! 482 ! first ABL level 526 !------------- 527 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 528 ! ! 6 *** MPI exchanges 529 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 530 ! 531 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1._wp, v_abl(:,:,:,nt_a) , 'T', -1._wp ) 532 CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1._wp , tq_abl(:,:,:,nt_a,jp_qa), 'T', 1._wp , kfillmode = jpfillnothing ) ! ++++ this should not be needed... 533 ! 534 #if defined key_iomput 535 ! 2D & first ABL level 536 IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh", pblh(:,: ) ) 483 537 IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl", u_abl(:,:,2,nt_a ) ) 484 538 IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl", v_abl(:,:,2,nt_a ) ) … … 489 543 IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta", pt_dta(:,:,2 ) ) 490 544 IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta", pq_dta(:,:,2 ) ) 491 ! all ABL levels 492 IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) ) 493 IF ( iom_use("v_abl" ) ) CALL iom_put ( "v_abl" , v_abl(:,:,2:jpka,nt_a ) ) 494 IF ( iom_use("t_abl" ) ) CALL iom_put ( "t_abl" , tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 495 IF ( iom_use("q_abl" ) ) CALL iom_put ( "q_abl" , tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 496 IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a ) ) 497 IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka ) ) 498 IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka ) ) 499 IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxl_abl(:,:,2:jpka ) ) 500 IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh" , pblh(:,: ) ) 501 ! debug (to be removed) 545 ! debug 2D 546 IF( ln_geos_winds ) THEN 547 IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) ) 548 IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) ) 549 END IF 550 IF( ln_hpgls_frc ) THEN 551 IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 552 IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 553 END IF 554 ! 3D (all ABL levels) 555 IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) ) 556 IF ( iom_use("v_abl" ) ) CALL iom_put ( "v_abl" , v_abl(:,:,2:jpka,nt_a ) ) 557 IF ( iom_use("t_abl" ) ) CALL iom_put ( "t_abl" , tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 558 IF ( iom_use("q_abl" ) ) CALL iom_put ( "q_abl" , tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 559 IF ( iom_use("tke_abl" ) ) CALL iom_put ( "tke_abl" , tke_abl(:,:,2:jpka,nt_a ) ) 560 IF ( iom_use("avm_abl" ) ) CALL iom_put ( "avm_abl" , avm_abl(:,:,2:jpka ) ) 561 IF ( iom_use("avt_abl" ) ) CALL iom_put ( "avt_abl" , avt_abl(:,:,2:jpka ) ) 562 IF ( iom_use("mxlm_abl") ) CALL iom_put ( "mxlm_abl", mxlm_abl(:,:,2:jpka ) ) 563 IF ( iom_use("mxld_abl") ) CALL iom_put ( "mxld_abl", mxld_abl(:,:,2:jpka ) ) 564 ! debug 3D 502 565 IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta", pu_dta(:,:,2:jpka) ) 503 566 IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta", pv_dta(:,:,2:jpka) ) 504 567 IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta", pt_dta(:,:,2:jpka) ) 505 568 IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta", pq_dta(:,:,2:jpka) ) 506 IF ( iom_use("coeft") ) CALL iom_put ( "coeft", z_cft(:,:,2:jpka) )507 569 IF( ln_geos_winds ) THEN 508 IF ( iom_use("u z1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) )509 IF ( iom_use("v z1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) )570 IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka) ) 571 IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgv_dta(:,:,2:jpka) ) 510 572 END IF 511 573 IF( ln_hpgls_frc ) THEN 512 IF ( iom_use("u z1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp))513 IF ( iom_use("v z1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp))574 IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) ) 575 IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", -pgv_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) ) 514 576 END IF 515 ! 577 #endif 516 578 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 517 579 ! ! 7 *** Finalize flux computation 518 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 519 580 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 581 ! 520 582 DO_2D_11_11 521 ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta ) 522 zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa ) 523 !zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ... 524 ! & ( R_dry*ztemp * ( 1._wp + rctv0*zhumi ) ) 525 zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 526 psen ( ji, jj ) = cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) 527 pevp ( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) ) 528 rhoa( ji, jj ) = zcff 583 ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta ) 584 zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa ) 585 zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 586 psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) !GS: negative sign to respect aerobulk convention 587 pevp( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) ) 588 rhoa( ji, jj ) = zcff 529 589 END_2D 530 590 531 591 DO_2D_01_01 532 zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac *( pssu(ji ,jj) + pssu(ji-1,jj) )533 zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac *( pssv(ji,jj ) + pssv(ji,jj-1) )592 zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * ( pssu(ji ,jj) + pssu(ji-1,jj) ) 593 zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * ( pssv(ji,jj ) + pssv(ji,jj-1) ) 534 594 END_2D 535 ! 595 ! 536 596 CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1.0_wp, zwnd_j(:,:) , 'T', -1.0_wp ) 537 597 ! … … 539 599 DO_2D_11_11 540 600 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 541 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) )! * msk_abl(ji,jj)601 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 542 602 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 543 603 544 604 pwndm (ji,jj) = zcff 545 605 ptaum (ji,jj) = zztmp * zcff … … 564 624 565 625 IF(sn_cfctl%l_prtctl) THEN 566 CALL prt_ctl( tab2d_1=p wndm , clinfo1=' abl_stp: wndm : ' )567 CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ')568 CALL prt_ctl( tab2d_ 2=ptauj , clinfo2= 'vtau: ' )626 CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ', mask1=umask, & 627 & tab2d_2=ptauj , clinfo2=' vtau : ', mask2=vmask ) 628 CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' ) 569 629 ENDIF 570 630 571 631 #if defined key_si3 572 ! ------------------------------------------------------------ ! 573 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 574 ! ------------------------------------------------------------ ! 575 DO_2D_00_00 576 577 zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 578 zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 579 580 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 581 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 582 & * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 583 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 584 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 585 & * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 586 END_2D 587 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp ) 588 ! 589 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' & 590 & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' ) 632 ! ------------------------------------------------------------ ! 633 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 634 ! ------------------------------------------------------------ ! 635 DO_2D_00_00 636 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) & 637 & * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) ) 638 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) & 639 & * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) ) 640 END_2D 641 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp ) 642 ! 643 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' & 644 & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' ) 645 ! ------------------------------------------------------------ ! 646 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 647 ! ------------------------------------------------------------ ! 648 DO_2D_00_00 649 650 zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 651 zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 652 653 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 654 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 655 & * ( zztmp1 - pssu_ice(ji,jj) ) 656 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 657 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 658 & * ( zztmp2 - pssv_ice(ji,jj) ) 659 END_2D 660 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice,'V', -1.0_wp ) 661 ! 662 IF(sn_cfctl%l_prtctl) THEN 663 CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask, & 664 & tab2d_2=ptauj_ice , clinfo2=' vtau_ice : ', mask2=vmask ) 665 END IF 591 666 #endif 592 667 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 599 674 END SUBROUTINE abl_stp 600 675 !=================================================================================================== 601 602 603 604 605 606 607 608 609 610 611 612 613 614 676 615 677 … … 634 696 !! (= Kz dz[Ub] * dz[Un] ) 635 697 !! --------------------------------------------------------------------- 636 INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn 698 INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn 637 699 INTEGER, DIMENSION(1:jpi ) :: ikbl 638 700 REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv 639 REAL(wp) :: zdU , zdV, zcff1,zshear,zbuoy,zsig, zustar2640 REAL(wp) :: zdU2, zdV2641 REAL(wp) :: zwndi, zwndj701 REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2 702 REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 ! zbuoy for BL89 703 REAL(wp) :: zwndi, zwndj 642 704 REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2 643 705 REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2 644 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF 706 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF 645 707 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a 646 708 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b … … 648 710 LOGICAL :: ln_Patankar = .FALSE. 649 711 LOGICAL :: ln_dumpvar = .FALSE. 650 LOGICAL , DIMENSION(1:jpi ) :: ln_foundl 712 LOGICAL , DIMENSION(1:jpi ) :: ln_foundl 651 713 ! 652 714 tind = nt_n … … 660 722 !------------- 661 723 ! 662 ! Compute vertical shear 724 ! Compute vertical shear 663 725 DO jk = 2, jpkam1 664 DO ji = 1, jpi665 zcff = 1.0_wp / e3w_abl( jk )**2 666 zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 726 DO ji = 1, jpi 727 zcff = 1.0_wp / e3w_abl( jk )**2 728 zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 667 729 zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 668 zsh2(ji,jk) = zdU+zdV 669 END DO 670 END DO 730 zsh2(ji,jk) = zdU+zdV !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 ) 731 END DO 732 END DO 671 733 ! 672 734 ! Compute brunt-vaisala frequency 673 735 DO jk = 2, jpkam1 674 DO ji = 1,jpi 675 zcff = grav * itvref / e3w_abl( jk ) 736 DO ji = 1,jpi 737 zcff = grav * itvref / e3w_abl( jk ) 676 738 zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta) 677 739 zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) & … … 679 741 zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi) 680 742 END DO 681 END DO 743 END DO 682 744 ! 683 745 ! Terms for the tridiagonal problem 684 746 DO jk = 2, jpkam1 685 DO ji = 1, jpi686 zshear = zsh2( ji, jk )! zsh2 is already multiplied by Avm_abl at this point687 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation688 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )689 690 z_elem_a( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal691 z_elem_c( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal747 DO ji = 1, jpi 748 zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point 749 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 750 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 751 752 z_elem_a( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk ) + Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal 753 z_elem_c( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk ) + Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal 692 754 IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE 693 z_elem_b( ji, jk )= e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &694 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) ! diagonal695 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )! right-hand-side755 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 756 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) ! diagonal 757 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) ) ! right-hand-side 696 758 ELSE 697 z_elem_b( ji, jk )= e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &698 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) & ! diagonal699 & - e3w_abl(jk) * rDt_abl * zbuoy700 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side759 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 760 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) & ! diagonal 761 & - e3w_abl(jk) * rDt_abl * zbuoy 762 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side 701 763 END IF 702 764 END DO 703 END DO 704 705 DO ji = 1,jpi ! vector opt. 706 zesrf = MAX( 4.63_wp * ustar2(ji,jj), tke_min ) 707 zetop = tke_min 708 z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp 709 z_elem_a ( ji, jpka ) = 0._wp; z_elem_c ( ji, jpka ) = 0._wp; z_elem_b ( ji, jpka ) = 1._wp 710 tke_abl( ji, jj, 1, nt_a ) = zesrf 711 tke_abl( ji, jj, jpka, nt_a ) = zetop 712 zbn2(ji,jj, 1) = zbn2( ji,jj, 2) 713 zsh2(ji, 1) = zsh2( ji, 2) 714 zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 715 zsh2(ji, jpka) = zsh2( ji , jpkam1) 716 END DO 765 END DO 766 767 DO ji = 1,jpi ! vector opt. 768 zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min ) 769 zetop = tke_min 770 771 z_elem_a ( ji, 1 ) = 0._wp 772 z_elem_c ( ji, 1 ) = 0._wp 773 z_elem_b ( ji, 1 ) = 1._wp 774 tke_abl ( ji, jj, 1, nt_a ) = zesrf 775 776 !++ Top Neumann B.C. 777 !z_elem_a ( ji, jpka ) = - 0.5 * rDt_abl * rn_Sch * (Avm_abl(ji,jj, jpka-1 )+Avm_abl(ji,jj, jpka )) / e3t_abl( jpka ) 778 !z_elem_c ( ji, jpka ) = 0._wp 779 !z_elem_b ( ji, jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 780 !tke_abl ( ji, jj, jpka, nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 781 782 !++ Top Dirichlet B.C. 783 z_elem_a ( ji, jpka ) = 0._wp 784 z_elem_c ( ji, jpka ) = 0._wp 785 z_elem_b ( ji, jpka ) = 1._wp 786 tke_abl ( ji, jj, jpka, nt_a ) = zetop 787 788 zbn2 ( ji, jj, 1 ) = zbn2 ( ji, jj, 2 ) 789 zsh2 ( ji, 1 ) = zsh2 ( ji, 2 ) 790 zbn2 ( ji, jj, jpka ) = zbn2 ( ji, jj, jpkam1 ) 791 zsh2 ( ji, jpka ) = zsh2 ( ji , jpkam1 ) 792 END DO 717 793 !! 718 794 !! Matrix inversion … … 720 796 DO ji = 1,jpi 721 797 zcff = 1._wp / z_elem_b( ji, 1 ) 722 zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 ) 723 tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a ) 724 END DO 725 726 DO jk = 2, jpka 798 zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 ) 799 tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a ) 800 END DO 801 802 DO jk = 2, jpka 727 803 DO ji = 1,jpi 728 804 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) ) … … 732 808 END DO 733 809 END DO 734 735 DO jk = jpkam1,1,-1 810 811 DO jk = jpkam1,1,-1 736 812 DO ji = 1,jpi 737 813 tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a) 738 814 END DO 739 815 END DO 740 741 !!FL should not be needed because of Patankar procedure 816 817 !!FL should not be needed because of Patankar procedure 742 818 tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min ) 743 819 … … 745 821 !! Diagnose PBL height 746 822 !! ---------------------------------------------------------- 747 748 749 ! 823 824 825 ! 750 826 ! arrays zRH, zFC and zCF are available at this point 751 827 ! and zFC(:, 1 ) = 0. 752 828 ! diagnose PBL height based on zsh2 and zbn2 753 829 zFC ( : ,1) = 0._wp 754 ikbl( 1:jpi ) = 0 755 830 ikbl( 1:jpi ) = 0 831 756 832 DO jk = 2,jpka 757 DO ji = 1, jpi 833 DO ji = 1, jpi 758 834 zcff = ghw_abl( jk-1 ) 759 835 zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) ) … … 781 857 ELSE 782 858 pblh( ji, jj ) = ghw_abl(jpka) 783 END IF 784 END DO 785 !------------- 786 END DO 787 !------------- 788 ! 789 ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 859 END IF 860 END DO 861 !------------- 862 END DO 863 !------------- 864 ! 865 ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 790 866 IF(ln_smth_pblh) THEN 791 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) 867 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing) 792 868 CALL smooth_pblh( pblh, msk_abl ) 793 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) 869 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing) 794 870 ENDIF 795 871 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 799 875 SELECT CASE ( nn_amxl ) 800 876 ! 801 CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom 802 # define zlup zRH 803 # define zldw zFC 877 CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom 878 # define zlup zRH 879 # define zldw zFC 804 880 DO jj = 1, jpj ! outer loop 805 881 ! 806 882 DO ji = 1, jpi 807 mxl_abl ( ji, jj, 1 ) = 0._wp 808 mxl_abl ( ji, jj, jpka ) = mxl_min 809 zldw( ji, 1 ) = 0._wp 810 zlup( ji, jpka ) = 0._wp 811 END DO 812 ! 813 DO jk = 2, jpkam1 814 DO ji = 1, jpi 815 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 816 mxl_abl( ji, jj, jk ) = MAX( mxl_min, & 817 & SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) ) 818 END DO 819 END DO 883 mxld_abl( ji, jj, 1 ) = mxl_min 884 mxld_abl( ji, jj, jpka ) = mxl_min 885 mxlm_abl( ji, jj, 1 ) = mxl_min 886 mxlm_abl( ji, jj, jpka ) = mxl_min 887 zldw ( ji, 1 ) = zrough(ji,jj) * rn_Lsfc 888 zlup ( ji, jpka ) = mxl_min 889 END DO 890 ! 891 DO jk = 2, jpkam1 892 DO ji = 1, jpi 893 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 894 mxlm_abl( ji, jj, jk ) = MAX( mxl_min, & 895 & SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) ) 896 END DO 897 END DO 820 898 ! 821 899 ! Limit mxl 822 DO jk = jpkam1,1,-1 823 DO ji = 1, jpi 824 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl _abl(ji, jj, jk) )825 END DO 826 END DO 900 DO jk = jpkam1,1,-1 901 DO ji = 1, jpi 902 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 903 END DO 904 END DO 827 905 ! 828 906 DO jk = 2, jpka 829 DO ji = 1, jpi 830 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) ) 831 END DO 832 END DO 907 DO ji = 1, jpi 908 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 909 END DO 910 END DO 911 ! 912 ! DO jk = 1, jpka 913 ! DO ji = 1, jpi 914 ! mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 915 ! mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ), zlup( ji, jk ) ) 916 ! END DO 917 ! END DO 833 918 ! 834 919 DO jk = 1, jpka 835 920 DO ji = 1, jpi 836 mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 837 END DO 838 END DO 839 ! 840 END DO 841 # undef zlup 842 # undef zldw 843 ! 844 ! 845 CASE ( 1 ) ! length-scale computed as the distance to the PBL height 846 DO jj = 1,jpj ! outer loop 847 ! 848 DO ji = 1, jpi ! vector opt. 849 zcff = 1._wp / pblh( ji, jj ) ! inverse of hbl 850 DO jk = 1, jpka 851 zsig = MIN( zcff * ghw_abl( jk ), 1.0_wp ) 852 zcff1 = pblh( ji, jj ) 853 mxl_abl( ji, jj, jk ) = mxl_min & 854 & + zsig * ( amx1*zcff1 + bmx1*mxl_min ) & 855 & + zsig * zsig * ( amx2*zcff1 + bmx2*mxl_min ) & 856 & + zsig**3 * ( amx3*zcff1 + bmx3*mxl_min ) & 857 & + zsig**4 * ( amx4*zcff1 + bmx4*mxl_min ) & 858 & + zsig**5 * ( amx5*zcff1 + bmx5*mxl_min ) 859 END DO 860 END DO 861 ! 862 END DO 921 ! zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 922 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 923 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 924 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 925 END DO 926 END DO 927 ! 928 END DO 929 # undef zlup 930 # undef zldw 931 ! 932 ! 933 CASE ( 1 ) ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom 934 # define zlup zRH 935 # define zldw zFC 936 DO jj = 1, jpj ! outer loop 937 ! 938 DO jk = 2, jpkam1 939 DO ji = 1,jpi 940 zcff = 1.0_wp / e3w_abl( jk )**2 941 zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 942 zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 943 zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 944 ENDDO 945 ENDDO 946 ! 947 DO ji = 1, jpi 948 zcff = zrough(ji,jj) * rn_Lsfc 949 mxld_abl ( ji, jj, 1 ) = zcff 950 mxld_abl ( ji, jj, jpka ) = mxl_min 951 mxlm_abl ( ji, jj, 1 ) = zcff 952 mxlm_abl ( ji, jj, jpka ) = mxl_min 953 zldw ( ji, 1 ) = zcff 954 zlup ( ji, jpka ) = mxl_min 955 END DO 956 ! 957 DO jk = 2, jpkam1 958 DO ji = 1, jpi 959 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 960 zcff = 2.0_wp*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) & 961 & + SQRT(rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.0_wp*zbuoy ) ) 962 mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff ) 963 END DO 964 END DO 965 ! 966 ! Limit mxl 967 DO jk = jpkam1,1,-1 968 DO ji = 1, jpi 969 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 970 END DO 971 END DO 972 ! 973 DO jk = 2, jpka 974 DO ji = 1, jpi 975 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 976 END DO 977 END DO 978 ! 979 DO jk = 1, jpka 980 DO ji = 1, jpi 981 !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 982 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 983 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 984 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 985 !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) ) 986 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 987 END DO 988 END DO 989 ! 990 END DO 991 # undef zlup 992 # undef zldw 863 993 ! 864 994 CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale 865 995 ! 866 # define zlup zRH 867 # define zldw zFC 996 # define zlup zRH 997 # define zldw zFC 868 998 ! zCF is used for matrix inversion 869 ! 999 ! 870 1000 DO jj = 1, jpj ! outer loop 871 872 DO ji = 1, jpi 873 zlup( ji, 1 ) = mxl_min 874 zldw( ji, 1 ) = mxl_min 1001 1002 DO ji = 1, jpi 1003 zcff = zrough(ji,jj) * rn_Lsfc 1004 zlup( ji, 1 ) = zcff 1005 zldw( ji, 1 ) = zcff 875 1006 zlup( ji, jpka ) = mxl_min 876 zldw( ji, jpka ) = mxl_min 877 END DO 878 1007 zldw( ji, jpka ) = mxl_min 1008 END DO 1009 879 1010 DO jk = 2,jpka-1 880 1011 DO ji = 1, jpi 881 1012 zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 882 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 883 END DO 884 END DO 1013 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 1014 END DO 1015 END DO 885 1016 !! 886 1017 !! BL89 search for lup 887 !! ---------------------------------------------------------- 888 DO jk=2,jpka-1 1018 !! ---------------------------------------------------------- 1019 DO jk=2,jpka-1 889 1020 ! 890 1021 DO ji = 1, jpi … … 892 1023 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 893 1024 ln_foundl(ji ) = .false. 894 END DO 895 ! 1025 END DO 1026 ! 896 1027 DO jkup=jk+1,jpka-1 897 1028 DO ji = 1, jpi 1029 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1030 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 898 1031 zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 899 & ( zb n2(ji,jj,jkup )*(ghw_abl(jkup )-ghw_abl(jk)) &900 & + zb n2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )1032 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & 1033 & + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) 901 1034 IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 902 1035 zcff2 = ghw_abl(jkup ) - ghw_abl(jk) 903 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 1036 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 904 1037 zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / & 905 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 906 zlup(ji,jk) = zcff 1038 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 1039 zlup(ji,jk) = zcff 1040 zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk) 907 1041 ln_foundl(ji) = .true. 908 1042 END IF … … 910 1044 END DO 911 1045 ! 912 END DO 1046 END DO 913 1047 !! 914 1048 !! BL89 search for ldwn 915 !! ---------------------------------------------------------- 916 DO jk=2,jpka-1 1049 !! ---------------------------------------------------------- 1050 DO jk=2,jpka-1 917 1051 ! 918 1052 DO ji = 1, jpi … … 920 1054 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 921 1055 ln_foundl(ji ) = .false. 922 END DO 923 ! 1056 END DO 1057 ! 924 1058 DO jkdwn=jk-1,1,-1 925 DO ji = 1, jpi 1059 DO ji = 1, jpi 1060 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1061 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 926 1062 zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) & 927 & * ( zb n2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &928 + zb n2(ji,jj,jkdwn )*(ghw_abl(jk)-ghw_abl(jkdwn )) )929 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1063 & * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 1064 + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) 1065 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 930 1066 zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 931 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 1067 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 932 1068 zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / & 933 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 934 zldw(ji,jk) = zcff 935 ln_foundl(ji) = .true. 936 END IF 937 END DO 938 END DO 939 ! 1069 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1070 zldw(ji,jk) = zcff 1071 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1072 ln_foundl(ji) = .true. 1073 END IF 1074 END DO 1075 END DO 1076 ! 940 1077 END DO 941 1078 942 1079 DO jk = 1, jpka 943 DO ji = 1, jpi 944 mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min ) 945 END DO 946 END DO 1080 DO ji = 1, jpi 1081 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1082 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1083 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1084 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1085 END DO 1086 END DO 947 1087 948 1088 END DO 949 # undef zlup 950 # undef zldw 951 ! 952 END SELECT 1089 # undef zlup 1090 # undef zldw 1091 ! 1092 CASE ( 3 ) ! Bougeault & Lacarrere 89 length-scale 1093 ! 1094 # define zlup zRH 1095 # define zldw zFC 1096 ! zCF is used for matrix inversion 1097 ! 1098 DO jj = 1, jpj ! outer loop 1099 ! 1100 DO jk = 2, jpkam1 1101 DO ji = 1,jpi 1102 zcff = 1.0_wp / e3w_abl( jk )**2 1103 zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 1104 zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 1105 zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 1106 ENDDO 1107 ENDDO 1108 zsh2(:, 1) = zsh2( :, 2) 1109 zsh2(:, jpka) = zsh2( :, jpkam1) 1110 1111 DO ji = 1, jpi 1112 zcff = zrough(ji,jj) * rn_Lsfc 1113 zlup( ji, 1 ) = zcff 1114 zldw( ji, 1 ) = zcff 1115 zlup( ji, jpka ) = mxl_min 1116 zldw( ji, jpka ) = mxl_min 1117 END DO 1118 1119 DO jk = 2,jpka-1 1120 DO ji = 1, jpi 1121 zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 1122 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 1123 END DO 1124 END DO 1125 !! 1126 !! BL89 search for lup 1127 !! ---------------------------------------------------------- 1128 DO jk=2,jpka-1 1129 ! 1130 DO ji = 1, jpi 1131 zCF(ji,1:jpka) = 0._wp 1132 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 1133 ln_foundl(ji ) = .false. 1134 END DO 1135 ! 1136 DO jkup=jk+1,jpka-1 1137 DO ji = 1, jpi 1138 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1139 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 1140 zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 1141 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & 1142 & + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) & 1143 & + 0.5_wp * e3t_abl(jkup) * rn_Rod * & 1144 & ( SQRT(tke_abl( ji, jj, jkup , nt_a ))*zsh2(ji,jkup ) & 1145 & + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) ) 1146 1147 IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1148 zcff2 = ghw_abl(jkup ) - ghw_abl(jk) 1149 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 1150 zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / & 1151 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 1152 zlup(ji,jk) = zcff 1153 zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk) 1154 ln_foundl(ji) = .true. 1155 END IF 1156 END DO 1157 END DO 1158 ! 1159 END DO 1160 !! 1161 !! BL89 search for ldwn 1162 !! ---------------------------------------------------------- 1163 DO jk=2,jpka-1 1164 ! 1165 DO ji = 1, jpi 1166 zCF(ji,1:jpka) = 0._wp 1167 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 1168 ln_foundl(ji ) = .false. 1169 END DO 1170 ! 1171 DO jkdwn=jk-1,1,-1 1172 DO ji = 1, jpi 1173 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1174 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 1175 zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) & 1176 & * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 1177 & +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) & 1178 & + 0.5_wp * e3t_abl(jkup) * rn_Rod * & 1179 & ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) & 1180 & + SQRT(tke_abl( ji, jj, jkdwn , nt_a ))*zsh2(ji,jkdwn ) ) 1181 1182 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1183 zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 1184 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 1185 zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / & 1186 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1187 zldw(ji,jk) = zcff 1188 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1189 ln_foundl(ji) = .true. 1190 END IF 1191 END DO 1192 END DO 1193 ! 1194 END DO 1195 1196 DO jk = 1, jpka 1197 DO ji = 1, jpi 1198 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1199 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1200 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1201 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1202 END DO 1203 END DO 1204 1205 END DO 1206 # undef zlup 1207 # undef zldw 1208 ! 1209 ! 1210 END SELECT 953 1211 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 954 1212 ! ! Finalize the computation of turbulent visc./diff. 955 1213 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 956 1214 957 1215 !------------- 958 1216 DO jj = 1, jpj ! outer loop 959 1217 !------------- 960 DO jk = 1, jpka 1218 DO jk = 1, jpka 961 1219 DO ji = 1, jpi ! vector opt. 962 zcff = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk ) &963 & * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) )964 zcff2 965 zcff = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )966 !!FL: MAX function probably useless because of the definition of mxl_min 1220 zcff = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk ) & 1221 & * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) ) 1222 zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z) 1223 zcff = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) ) 1224 !!FL: MAX function probably useless because of the definition of mxl_min 967 1225 Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak ) 968 Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak ) 969 END DO 970 END DO 971 !------------- 972 END DO 1226 Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak ) 1227 END DO 1228 END DO 1229 !------------- 1230 END DO 973 1231 !------------- 974 1232 … … 988 1246 !! 989 1247 !! --------------------------------------------------------------------- 990 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk 991 1248 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk 1249 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d 992 1250 INTEGER :: ji,jj 993 994 995 1251 REAL(wp) :: smth_a, smth_b 1252 REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY 1253 REAL(wp) :: zumsk,zvmsk 996 1254 !! 997 1255 !!========================================================= … … 1005 1263 zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk 1006 1264 END_2D 1007 1008 1265 1266 DO_2D_10_11 1009 1267 zvmsk = msk(ji,jj) * msk(ji,jj+1) 1010 1268 zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk 1011 1012 1013 1269 END_2D 1270 1271 DO_2D_10_00 1014 1272 zFY ( ji, jj ) = zdY ( ji, jj ) & 1015 1273 & + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) & 1016 1274 & - (zdX ( ji, jj ) - zdX( ji-1, jj )) ) 1017 1275 END_2D 1018 1276 1019 1277 DO_2D_00_10 … … 1029 1287 & +zFY( ji, jj ) - zFY( ji, jj-1 ) ) 1030 1288 END_2D 1031 !! 1289 1032 1290 !--------------------------------------------------------------------------------------------------- 1033 1291 END SUBROUTINE smooth_pblh -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ABL/ablrst.F90
r13135 r13220 109 109 CALL iom_delay_rst( 'WRITE', 'ABL', numraw ) ! save only abl delayed global communication variables 110 110 111 ! Prognostic variables111 ! Prognostic (after timestep + swap time indices = now timestep) variables 112 112 CALL iom_rstput( iter, nitrst, numraw, 'u_abl', u_abl(:,:,:,nt_n ) ) 113 113 CALL iom_rstput( iter, nitrst, numraw, 'v_abl', v_abl(:,:,:,nt_n ) ) … … 117 117 CALL iom_rstput( iter, nitrst, numraw, 'avm_abl', avm_abl(:,:,: ) ) 118 118 CALL iom_rstput( iter, nitrst, numraw, 'avt_abl', avt_abl(:,:,: ) ) 119 CALL iom_rstput( iter, nitrst, numraw, 'mxl_abl', mxl_abl(:,:,: ) )119 CALL iom_rstput( iter, nitrst, numraw,'mxld_abl',mxld_abl(:,:,: ) ) 120 120 CALL iom_rstput( iter, nitrst, numraw, 'pblh', pblh(:,: ) ) 121 121 ! … … 172 172 CALL iom_get( numrar, jpdom_autoglo, 'avm_abl', avm_abl(:,:,: ) ) 173 173 CALL iom_get( numrar, jpdom_autoglo, 'avt_abl', avt_abl(:,:,: ) ) 174 CALL iom_get( numrar, jpdom_autoglo, 'mxl_abl', mxl_abl(:,:,: ) )174 CALL iom_get( numrar, jpdom_autoglo,'mxld_abl',mxld_abl(:,:,: ) ) 175 175 CALL iom_get( numrar, jpdom_autoglo, 'pblh', pblh(:,: ) ) 176 176 CALL iom_delay_rst( 'READ', 'ABL', numrar ) ! read only abl delayed global communication variables -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ABL/par_abl.F90
r13135 r13220 28 28 LOGICAL , PUBLIC :: ln_hpgls_frc !: forcing of ABL winds by large-scale pressure gradient 29 29 LOGICAL , PUBLIC :: ln_smth_pblh !: smoothing of atmospheric PBL height 30 !LOGICAL , PUBLIC :: ln_topbc_neumann = .FALSE. !: idealised testcases only 30 31 31 LOGICAL , PUBLIC :: ln_rstart_abl !: (de)activate abl restart32 CHARACTER(len=256), PUBLIC :: cn_ablrst_in !: suffix of abl restart name (input)33 CHARACTER(len=256), PUBLIC :: cn_ablrst_out !: suffix of abl restart name (output)34 CHARACTER(len=256), PUBLIC :: cn_ablrst_indir !: abl restart input directory35 CHARACTER(len=256), PUBLIC :: cn_ablrst_outdir !: abl restart output directory32 LOGICAL , PUBLIC :: ln_rstart_abl !: (de)activate abl restart 33 CHARACTER(len=256), PUBLIC :: cn_ablrst_in !: suffix of abl restart name (input) 34 CHARACTER(len=256), PUBLIC :: cn_ablrst_out !: suffix of abl restart name (output) 35 CHARACTER(len=256), PUBLIC :: cn_ablrst_indir !: abl restart input directory 36 CHARACTER(len=256), PUBLIC :: cn_ablrst_outdir !: abl restart output directory 36 37 37 38 !!--------------------------------------------------------------------- … … 46 47 REAL(wp), PUBLIC, PARAMETER :: rn_Cek = 258._wp !: Ekman constant for Richardson number 47 48 REAL(wp), PUBLIC, PARAMETER :: rn_epssfc = 1._wp / ( 1._wp + 2.8_wp * 2.8_wp ) 48 REAL(wp), PUBLIC :: rn_ ceps !: namelist parameter49 REAL(wp), PUBLIC :: rn_ cm !: namelist parameter50 REAL(wp), PUBLIC :: rn_ ct !: namelist parameter51 REAL(wp), PUBLIC :: rn_ ce !: namelist parameter49 REAL(wp), PUBLIC :: rn_Ceps !: namelist parameter 50 REAL(wp), PUBLIC :: rn_Cm !: namelist parameter 51 REAL(wp), PUBLIC :: rn_Ct !: namelist parameter 52 REAL(wp), PUBLIC :: rn_Ce !: namelist parameter 52 53 REAL(wp), PUBLIC :: rn_Rod !: namelist parameter 53 54 REAL(wp), PUBLIC :: rn_Sch 55 REAL(wp), PUBLIC :: rn_Esfc 56 REAL(wp), PUBLIC :: rn_Lsfc 54 57 REAL(wp), PUBLIC :: mxl_min 55 58 REAL(wp), PUBLIC :: rn_ldyn_min !: namelist parameter -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ABL/sbcabl.F90
r13135 r13220 71 71 & ln_hpgls_frc, ln_geos_winds, nn_dyn_restore, & 72 72 & rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max, & 73 & nn_amxl, rn_ cm, rn_ct, rn_ce, rn_ceps, rn_Rod, rn_Ric, &73 & nn_amxl, rn_Cm, rn_Ct, rn_Ce, rn_Ceps, rn_Rod, rn_Ric, & 74 74 & ln_smth_pblh 75 75 !!--------------------------------------------------------------------- 76 76 77 ! Namelist namsbc_abl in reference namelist : ABL parameters77 ! Namelist namsbc_abl in reference namelist : ABL parameters 78 78 READ ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 ) 79 79 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist' ) 80 ! Namelist namsbc_abl in configuration namelist : ABL parameters 80 ! 81 ! Namelist namsbc_abl in configuration namelist : ABL parameters 81 82 READ ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 ) 82 83 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist' ) … … 165 166 rn_Sch = rn_ce / rn_cm 166 167 mxl_min = (avm_bak / rn_cm) / sqrt( tke_min ) 168 rn_Esfc = 1._wp / SQRT(rn_cm*rn_ceps) 169 rn_Lsfc = vkarmn * SQRT(SQRT(rn_cm*rn_ceps)) / rn_cm 167 170 168 171 IF(lwp) THEN … … 171 174 WRITE(numout,*) ' ~~~~~~~~~~~' 172 175 IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 173 IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height ' 176 IF(nn_amxl==1) WRITE(numout,*) 'Modified Deardorff 80 length-scale ' 177 IF(nn_amxl==2) WRITE(numout,*) 'Bougeault and Lacarrere length-scale ' 178 IF(nn_amxl==3) WRITE(numout,*) 'Rodier et al. length-scale ' 174 179 WRITE(numout,*) ' Minimum value of atmospheric TKE = ',tke_min,' m^2 s^-2' 175 180 WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m' … … 178 183 WRITE(numout,*) ' Constant for Schmidt number = ',rn_Sch 179 184 WRITE(numout,*) ' Constant for TKE dissipation = ',rn_Ceps 185 WRITE(numout,*) ' Constant for TKE sfc boundary condition = ',rn_Esfc 186 WRITE(numout,*) ' Constant for mxl sfc boundary condition = ',rn_Lsfc 180 187 END IF 181 188 … … 202 209 ! ABL timestep 203 210 rDt_abl = nn_fsbc * rn_Dt 211 IF(lwp) WRITE(numout,*) ' ABL timestep = ', rDt_abl,' s' 204 212 205 213 ! Check parameters for dynamics … … 248 256 zcff = 2._wp * omega * SIN( rad * 90._wp ) !++ fmax 249 257 rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 250 !!GS: alternative shape251 !rest_eq(:,:) = SIN( 0.5_wp*rpi*(zcff - ABS(ff_t(:,:))) / (zcff - 3.e-5) )**8252 !WHERE(ABS(ff_t(:,:)).LE.3.e-5) rest_eq(:,:) = 1._wp253 258 ELSE 254 259 rest_eq(:,:) = 1._wp … … 271 276 CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 272 277 273 u_abl(:,:,:,nt_n ) = sf(jp_wndi)%fnow(:,:,:)274 v_abl(:,:,:,nt_n ) = sf(jp_wndj)%fnow(:,:,:)278 u_abl(:,:,:,nt_n ) = sf(jp_wndi)%fnow(:,:,:) 279 v_abl(:,:,:,nt_n ) = sf(jp_wndj)%fnow(:,:,:) 275 280 tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 276 281 tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) … … 279 284 avm_abl(:,:,: ) = avm_bak 280 285 avt_abl(:,:,: ) = avt_bak 281 mxl_abl(:,:,: ) = mxl_min282 286 pblh (:,: ) = ghw_abl( 3 ) !<-- assume that the pbl contains 3 grid points 283 287 u_abl (:,:,:,nt_a ) = 0._wp … … 285 289 tq_abl (:,:,:,nt_a,: ) = 0._wp 286 290 tke_abl(:,:,:,nt_a ) = 0._wp 291 292 mxlm_abl(:,:,: ) = mxl_min 293 mxld_abl(:,:,: ) = mxl_min 287 294 ENDIF 288 295 … … 335 342 & tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in 336 343 & sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in 344 & sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in 337 345 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in 338 346 & tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/iceistate.F90
r13135 r13220 32 32 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 33 33 USE fldread ! read input fields 34 35 # if defined key_agrif 36 USE agrif_oce 37 USE agrif_ice 38 USE agrif_ice_interp 39 # endif 34 40 35 41 IMPLICIT NONE … … 168 174 ! 2) overwrite some of the fields with namelist parameters or netcdf file 169 175 !------------------------------------------------------------------------ 176 177 170 178 IF( ln_iceini ) THEN 171 179 ! !---------------! 172 IF( ln_iceini_file )THEN ! Read a file ! 173 ! !---------------! 174 WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp 175 ELSEWHERE ; zswitch(:,:) = 0._wp 180 181 IF( Agrif_Root() ) THEN 182 183 IF( ln_iceini_file )THEN ! Read a file ! 184 ! !---------------! 185 WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp 186 ELSEWHERE ; zswitch(:,:) = 0._wp 187 END WHERE 188 ! 189 CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step 190 ! 191 ! -- mandatory fields -- ! 192 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1) 193 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1) 194 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1) 195 196 ! -- optional fields -- ! 197 ! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature) 198 ! 199 ! ice salinity 200 IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) & 201 & si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 202 ! 203 ! temperatures 204 IF ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. & 205 & TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN 206 si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 207 si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 208 si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 209 ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 210 si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 211 ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 212 si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 213 ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s 214 si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 215 ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i 216 si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 217 ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su 218 si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 219 ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i 220 si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 221 ENDIF 222 ! 223 ! pond concentration 224 IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 225 & si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 226 & * si(jp_ati)%fnow(:,:,1) 227 ! 228 ! pond depth 229 IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 230 & si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 231 ! 232 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1) 233 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1) 234 zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1) 235 ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1) 236 zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1) 237 zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1) 238 ! 239 ! change the switch for the following 240 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 241 ELSEWHERE ; zswitch(:,:) = 0._wp 242 END WHERE 243 244 ! !---------------! 245 ELSE ! Read namelist ! 246 ! !---------------! 247 ! no ice if (sst - Tfreez) >= thresold 248 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 249 ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) 250 END WHERE 251 ! 252 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 253 WHERE( ff_t(:,:) >= 0._wp ) 254 zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) 255 zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 256 zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 257 zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 258 ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 259 zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 260 ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 261 zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 262 zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 263 ELSEWHERE 264 zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 265 zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 266 zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 267 zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 268 ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 269 zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 270 ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 271 zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 272 zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 273 END WHERE 274 ! 275 ENDIF 276 277 278 279 ! make sure ponds = 0 if no ponds scheme 280 IF ( .NOT.ln_pnd ) THEN 281 zapnd_ini(:,:) = 0._wp 282 zhpnd_ini(:,:) = 0._wp 283 ENDIF 284 285 !-------------! 286 ! fill fields ! 287 !-------------! 288 ! select ice covered grid points 289 npti = 0 ; nptidx(:) = 0 290 DO_2D_11_11 291 IF ( zht_i_ini(ji,jj) > 0._wp ) THEN 292 npti = npti + 1 293 nptidx(npti) = (jj - 1) * jpi + ji 294 ENDIF 295 END_2D 296 297 ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj) 298 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini ) 299 CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini ) 300 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini ) 301 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini ) 302 CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini ) 303 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini ) 304 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini ) 305 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini ) 306 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini ) 307 308 ! allocate temporary arrays 309 ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 310 & zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 311 312 ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 313 CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), & 314 & zhi_2d , zhs_2d , zai_2d , & 315 & t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), & 316 & zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d ) 317 318 ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 319 DO jl = 1, jpl 320 zti_3d(:,:,jl) = rt0 * tmask(:,:,1) 321 zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 322 END DO 323 CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i ) 324 CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s ) 325 CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i ) 326 CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d ) 327 CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d ) 328 CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su ) 329 CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i ) 330 CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip ) 331 CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip ) 332 333 ! deallocate temporary arrays 334 DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 335 & zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d ) 336 337 ! calculate extensive and intensive variables 338 CALL ice_var_salprof ! for sz_i 339 DO jl = 1, jpl 340 DO_2D_11_11 341 v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 342 v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 343 sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl) 344 END_2D 345 END DO 346 ! 347 DO jl = 1, jpl 348 DO_3D_11_11( 1, nlay_s ) 349 t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl) 350 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * & 351 & rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 352 END_3D 353 END DO 354 ! 355 DO jl = 1, jpl 356 DO_3D_11_11( 1, nlay_i ) 357 t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 358 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 359 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 360 & rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + & 361 & rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) & 362 & - rcp * ( ztmelts - rt0 ) ) 363 END_3D 364 END DO 365 366 ! Melt ponds 367 WHERE( a_i > epsi10 ) 368 a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 369 ELSEWHERE 370 a_ip_frac(:,:,:) = 0._wp 176 371 END WHERE 372 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 373 374 ! specific temperatures for coupled runs 375 tn_ice(:,:,:) = t_su(:,:,:) 376 t1_ice(:,:,:) = t_i (:,:,1,:) 177 377 ! 178 CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step 179 ! 180 ! -- mandatory fields -- ! 181 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1) 182 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1) 183 zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1) 184 185 ! -- optional fields -- ! 186 ! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature) 187 ! 188 ! ice salinity 189 IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) & 190 & si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 191 ! 192 ! temperatures 193 IF ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. & 194 & TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN 195 si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 196 si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 197 si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 198 ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 199 si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 200 ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 201 si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 202 ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s 203 si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 204 ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i 205 si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 206 ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su 207 si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 208 ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i 209 si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 210 ENDIF 211 ! 212 ! pond concentration 213 IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 214 & si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 215 & * si(jp_ati)%fnow(:,:,1) 216 ! 217 ! pond depth 218 IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 219 & si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 220 ! 221 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1) 222 ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1) 223 zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1) 224 ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1) 225 zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1) 226 zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1) 227 ! 228 ! change the switch for the following 229 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 230 ELSEWHERE ; zswitch(:,:) = 0._wp 378 379 #if defined key_agrif 380 ELSE 381 382 Agrif_SpecialValue = -9999. 383 Agrif_UseSpecialValue = .TRUE. 384 CALL Agrif_init_variable(tra_iceini_id,procname=interp_tra_ice) 385 use_sign_north = .TRUE. 386 sign_north = -1. 387 CALL Agrif_init_variable(u_iceini_id ,procname=interp_u_ice) 388 CALL Agrif_init_variable(v_iceini_id ,procname=interp_v_ice) 389 Agrif_SpecialValue = 0._wp 390 use_sign_north = .FALSE. 391 Agrif_UseSpecialValue = .FALSE. 392 ! lbc ???? 393 ! Here we know : a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, t_su, e_s, e_i 394 CALL ice_var_glo2eqv 395 CALL ice_var_zapsmall 396 CALL ice_var_agg(2) 397 398 ! Melt ponds 399 WHERE( a_i > epsi10 ) 400 a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 401 ELSEWHERE 402 a_ip_frac(:,:,:) = 0._wp 231 403 END WHERE 232 ! !---------------! 233 ELSE ! Read namelist ! 234 ! !---------------! 235 ! no ice if (sst - Tfreez) >= thresold 236 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 237 ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) 238 END WHERE 239 ! 240 ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 241 WHERE( ff_t(:,:) >= 0._wp ) 242 zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) 243 zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 244 zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 245 zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 246 ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 247 zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 248 ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 249 zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 250 zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 404 WHERE( a_ip > 0._wp ) ! ??????? 405 h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 251 406 ELSEWHERE 252 zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 253 zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 254 zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 255 zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 256 ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 257 zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 258 ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 259 zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 260 zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 261 END WHERE 262 ! 263 ENDIF 264 265 ! make sure ponds = 0 if no ponds scheme 266 IF ( .NOT.ln_pnd ) THEN 267 zapnd_ini(:,:) = 0._wp 268 zhpnd_ini(:,:) = 0._wp 269 ENDIF 270 271 !-------------! 272 ! fill fields ! 273 !-------------! 274 ! select ice covered grid points 275 npti = 0 ; nptidx(:) = 0 276 DO_2D_11_11 277 IF ( zht_i_ini(ji,jj) > 0._wp ) THEN 278 npti = npti + 1 279 nptidx(npti) = (jj - 1) * jpi + ji 280 ENDIF 281 END_2D 282 283 ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj) 284 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini ) 285 CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini ) 286 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini ) 287 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini ) 288 CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini ) 289 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini ) 290 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini ) 291 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini ) 292 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini ) 293 294 ! allocate temporary arrays 295 ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 296 & zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 297 298 ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 299 CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), & 300 & zhi_2d , zhs_2d , zai_2d , & 301 & t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), & 302 & zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d ) 303 304 ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 305 DO jl = 1, jpl 306 zti_3d(:,:,jl) = rt0 * tmask(:,:,1) 307 zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 308 END DO 309 CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i ) 310 CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s ) 311 CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i ) 312 CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d ) 313 CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d ) 314 CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su ) 315 CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i ) 316 CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip ) 317 CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip ) 318 319 ! deallocate temporary arrays 320 DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 321 & zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d ) 322 323 ! calculate extensive and intensive variables 324 CALL ice_var_salprof ! for sz_i 325 DO jl = 1, jpl 326 DO_2D_11_11 327 v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 328 v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 329 sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl) 330 END_2D 331 END DO 332 ! 333 DO jl = 1, jpl 334 DO_3D_11_11( 1, nlay_s ) 335 t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl) 336 e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * & 337 & rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 338 END_3D 339 END DO 340 ! 341 DO jl = 1, jpl 342 DO_3D_11_11( 1, nlay_i ) 343 t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 344 ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 345 e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 346 & rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + & 347 & rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) & 348 & - rcp * ( ztmelts - rt0 ) ) 349 END_3D 350 END DO 351 352 ! Melt ponds 353 WHERE( a_i > epsi10 ) 354 a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 355 ELSEWHERE 356 a_ip_frac(:,:,:) = 0._wp 357 END WHERE 358 v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 359 360 ! specific temperatures for coupled runs 361 tn_ice(:,:,:) = t_su(:,:,:) 362 t1_ice(:,:,:) = t_i (:,:,1,:) 363 ! 407 h_ip(:,:,:) = 0._wp 408 END WHERE 409 410 tn_ice(:,:,:) = t_su(:,:,:) 411 t1_ice(:,:,:) = t_i (:,:,1,:) 412 #endif 413 ENDIF ! Agrif_Root 364 414 ENDIF ! ln_iceini 365 415 ! -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icestp.F90
r12489 r13220 240 240 CALL par_init ! set some ice run parameters 241 241 ! 242 #if defined key_agrif 243 CALL Agrif_Declare_Var_ice ! " " " " " Sea ice 244 #endif 245 ! 242 246 ! ! Allocate the ice arrays (sbc_ice already allocated in sbc_init) 243 247 ierr = ice_alloc () ! ice variables -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/NST/agrif_ice.F90
r10068 r13220 16 16 17 17 INTEGER, PUBLIC :: u_ice_id, v_ice_id, tra_ice_id 18 INTEGER, PUBLIC :: u_iceini_id, v_iceini_id, tra_iceini_id 18 19 INTEGER, PUBLIC :: nbstep_ice = 0 ! child time position in sea-ice model 19 20 -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/NST/agrif_ice_interp.F90
r10069 r13220 14 14 !!---------------------------------------------------------------------- 15 15 !! agrif_interp_ice : interpolation of ice at "after" sea-ice time step 16 !! agrif_interp_u_ice : atomic routine to interpolate u_ice17 !! agrif_interp_v_ice : atomic routine to interpolate v_ice18 !! agrif_interp_tra_ice : atomic routine to interpolate ice properties16 !! interp_u_ice : atomic routine to interpolate u_ice 17 !! interp_v_ice : atomic routine to interpolate v_ice 18 !! interp_tra_ice : atomic routine to interpolate ice properties 19 19 !!---------------------------------------------------------------------- 20 20 USE par_oce … … 23 23 USE ice 24 24 USE agrif_ice 25 USE agrif_oce 25 26 USE phycst , ONLY: rt0 26 27 … … 29 30 30 31 PUBLIC agrif_interp_ice ! called by agrif_user.F90 32 PUBLIC interp_tra_ice, interp_u_ice, interp_v_ice ! called by iceistate.F90 31 33 32 34 !!---------------------------------------------------------------------- … … 68 70 Agrif_SpecialValue = -9999. 69 71 Agrif_UseSpecialValue = .TRUE. 72 73 use_sign_north = .TRUE. 74 sign_north = -1. 75 if (cd_type == 'T') use_sign_north = .FALSE. 76 70 77 SELECT CASE( cd_type ) 71 78 CASE('U') ; CALL Agrif_Bc_variable( u_ice_id , procname=interp_u_ice , calledweight=zbeta ) … … 75 82 Agrif_SpecialValue = 0._wp 76 83 Agrif_UseSpecialValue = .FALSE. 84 85 use_sign_north = .FALSE. 77 86 ! 78 87 END SUBROUTINE agrif_interp_ice … … 156 165 ! and it is ok since we conserve tracers (same as in the ocean). 157 166 ALLOCATE( ztab(SIZE(a_i,1),SIZE(a_i,2),SIZE(ptab,3)) ) 158 167 159 168 IF( before ) THEN ! parent grid 160 169 jm = 1 -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/NST/agrif_ice_update.F90
r12377 r13220 66 66 CALL Agrif_Update_Variable( tra_ice_id , locupdate=(/1,0/), procname = update_tra_ice ) 67 67 #endif 68 use_sign_north = .TRUE. 69 sign_north = -1. 70 68 71 # if ! defined DECAL_FEEDBACK 69 72 CALL Agrif_Update_Variable( u_ice_id , procname = update_u_ice ) … … 73 76 CALL Agrif_Update_Variable( v_ice_id , locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname=update_v_ice) 74 77 #endif 78 use_sign_north = .FALSE. 75 79 ! CALL Agrif_Update_Variable( tra_ice_id , locupdate=(/0,2/), procname = update_tra_ice ) 76 80 ! CALL Agrif_Update_Variable( u_ice_id , locupdate=(/0,1/), procname = update_u_ice ) -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/NST/agrif_oce.F90
r12377 r13220 19 19 20 20 ! !!* Namelist namagrif: AGRIF parameters 21 LOGICAL , PUBLIC :: ln_init_chfrpar = .FALSE. !: set child grids initial state from parent 21 22 LOGICAL , PUBLIC :: ln_agrif_2way = .TRUE. !: activate two way nesting 22 23 LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !: use zeros (.false.) or not (.true.) in … … 29 30 ! 30 31 INTEGER , PUBLIC, PARAMETER :: nn_sponge_len = 2 !: Sponge width (in number of parent grid points) 32 31 33 LOGICAL , PUBLIC :: spongedoneT = .FALSE. !: tracer sponge layer indicator 32 34 LOGICAL , PUBLIC :: spongedoneU = .FALSE. !: dynamics sponge layer indicator … … 49 51 INTEGER , PUBLIC, SAVE :: Kbb_a, Kmm_a, Krhs_a !: AGRIF module-specific copies of time-level indices 50 52 51 # if defined key_vertical52 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent 53 54 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt_parent, mbku_parent, mbkv_parent 54 # endif55 55 56 56 INTEGER, PUBLIC :: tsn_id ! AGRIF profile for tracers interpolation and update … … 58 58 INTEGER, PUBLIC :: un_update_id, vn_update_id ! AGRIF profiles for udpates 59 59 INTEGER, PUBLIC :: tsn_sponge_id, un_sponge_id, vn_sponge_id ! AGRIF profiles for sponge layers 60 INTEGER, PUBLIC :: tsini_id, uini_id, vini_id, sshini_id ! AGRIF profile for initialization 60 61 # if defined key_top 61 62 INTEGER, PUBLIC :: trn_id, trn_sponge_id … … 68 69 INTEGER, PUBLIC :: mbkt_id, ht0_id 69 70 INTEGER, PUBLIC :: kindic_agr 71 72 ! North fold 73 !$AGRIF_DO_NOT_TREAT 74 LOGICAL, PUBLIC :: use_sign_north 75 REAL, PUBLIC :: sign_north 76 LOGICAL, PUBLIC :: l_ini_child = .FALSE. 77 # if defined key_vertical 78 LOGICAL, PUBLIC :: l_vremap = .TRUE. 79 # else 80 LOGICAL, PUBLIC :: l_vremap = .FALSE. 81 # endif 82 !$AGRIF_END_DO_NOT_TREAT 70 83 71 84 !!---------------------------------------------------------------------- … … 91 104 & tabspongedone_trn(jpi,jpj), & 92 105 # endif 93 # if defined key_vertical94 106 & ht0_parent(jpi,jpj), mbkt_parent(jpi,jpj), & 95 107 & hu0_parent(jpi,jpj), mbku_parent(jpi,jpj), & 96 108 & hv0_parent(jpi,jpj), mbkv_parent(jpi,jpj), & 97 # endif98 109 & tabspongedone_u (jpi,jpj), & 99 110 & tabspongedone_v (jpi,jpj), STAT = ierr(1) ) -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/NST/agrif_oce_interp.F90
r12377 r13220 34 34 USE lib_mpp 35 35 USE vremap 36 USE lbclnk 36 37 37 38 IMPLICIT NONE … … 44 45 PUBLIC interpunb, interpvnb , interpub2b, interpvb2b 45 46 PUBLIC interpe3t 46 #if defined key_vertical47 47 PUBLIC interpht0, interpmbkt 48 # endif 48 PUBLIC agrif_initts, agrif_initssh 49 49 50 INTEGER :: bdy_tinterp = 0 50 51 … … 89 90 Agrif_UseSpecialValue = ln_spc_dyn 90 91 ! 92 use_sign_north = .TRUE. 93 sign_north = -1. 91 94 CALL Agrif_Bc_variable( un_interp_id, procname=interpun ) 92 95 CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn ) 96 use_sign_north = .FALSE. 93 97 ! 94 98 Agrif_UseSpecialValue = .FALSE. 95 99 ! 96 100 ! --- West --- ! 97 ibdy1 = 2 98 ibdy2 = 1+nbghostcells 99 ! 100 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 101 IF( lk_west ) THEN 102 ibdy1 = 2 103 ibdy2 = 1+nbghostcells 104 ! 105 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 106 DO ji = mi0(ibdy1), mi1(ibdy2) 107 uu_b(ji,:,Krhs_a) = 0._wp 108 109 DO jk = 1, jpkm1 110 DO jj = 1, jpj 111 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 112 END DO 113 END DO 114 115 DO jj = 1, jpj 116 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 117 END DO 118 END DO 119 ENDIF 120 ! 101 121 DO ji = mi0(ibdy1), mi1(ibdy2) 102 uu_b(ji,:,Krhs_a) = 0._wp 103 122 zub(ji,:) = 0._wp ! Correct transport 104 123 DO jk = 1, jpkm1 105 124 DO jj = 1, jpj 106 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 107 END DO 108 END DO 109 110 DO jj = 1, jpj 111 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 112 END DO 113 END DO 114 ENDIF 115 ! 116 DO ji = mi0(ibdy1), mi1(ibdy2) 117 zub(ji,:) = 0._wp ! Correct transport 118 DO jk = 1, jpkm1 119 DO jj = 1, jpj 120 zub(ji,jj) = zub(ji,jj) & 121 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a)*umask(ji,jj,jk) 122 END DO 123 END DO 124 DO jj=1,jpj 125 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 126 END DO 127 128 DO jk = 1, jpkm1 129 DO jj = 1, jpj 130 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a)-zub(ji,jj)) * umask(ji,jj,jk) 131 END DO 132 END DO 133 END DO 134 135 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 136 DO ji = mi0(ibdy1), mi1(ibdy2) 137 zvb(ji,:) = 0._wp 125 zub(ji,jj) = zub(ji,jj) & 126 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a)*umask(ji,jj,jk) 127 END DO 128 END DO 129 DO jj=1,jpj 130 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 131 END DO 132 138 133 DO jk = 1, jpkm1 139 134 DO jj = 1, jpj 140 zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 141 END DO 142 END DO 143 DO jj = 1, jpj 144 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 145 END DO 135 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a)-zub(ji,jj)) * umask(ji,jj,jk) 136 END DO 137 END DO 138 END DO 139 140 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 141 DO ji = mi0(ibdy1), mi1(ibdy2) 142 zvb(ji,:) = 0._wp 143 DO jk = 1, jpkm1 144 DO jj = 1, jpj 145 zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 146 END DO 147 END DO 148 DO jj = 1, jpj 149 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 150 END DO 151 DO jk = 1, jpkm1 152 DO jj = 1, jpj 153 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a)-zvb(ji,jj))*vmask(ji,jj,jk) 154 END DO 155 END DO 156 END DO 157 ENDIF 158 ENDIF 159 160 ! --- East --- ! 161 IF( lk_east) THEN 162 ibdy1 = jpiglo-1-nbghostcells 163 ibdy2 = jpiglo-2 164 ! 165 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 166 DO ji = mi0(ibdy1), mi1(ibdy2) 167 uu_b(ji,:,Krhs_a) = 0._wp 168 DO jk = 1, jpkm1 169 DO jj = 1, jpj 170 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) & 171 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 172 END DO 173 END DO 174 DO jj = 1, jpj 175 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 176 END DO 177 END DO 178 ENDIF 179 ! 180 DO ji = mi0(ibdy1), mi1(ibdy2) 181 zub(ji,:) = 0._wp ! Correct transport 146 182 DO jk = 1, jpkm1 147 183 DO jj = 1, jpj 148 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a)-zvb(ji,jj))*vmask(ji,jj,jk) 149 END DO 150 END DO 151 END DO 152 ENDIF 153 154 ! --- East --- ! 155 ibdy1 = jpiglo-1-nbghostcells 156 ibdy2 = jpiglo-2 157 ! 158 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 159 DO ji = mi0(ibdy1), mi1(ibdy2) 160 uu_b(ji,:,Krhs_a) = 0._wp 184 zub(ji,jj) = zub(ji,jj) & 185 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 186 END DO 187 END DO 188 DO jj=1,jpj 189 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 190 END DO 191 161 192 DO jk = 1, jpkm1 162 193 DO jj = 1, jpj 163 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) & 164 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 165 END DO 166 END DO 167 DO jj = 1, jpj 168 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 169 END DO 170 END DO 171 ENDIF 172 ! 173 DO ji = mi0(ibdy1), mi1(ibdy2) 174 zub(ji,:) = 0._wp ! Correct transport 175 DO jk = 1, jpkm1 176 DO jj = 1, jpj 177 zub(ji,jj) = zub(ji,jj) & 178 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 179 END DO 180 END DO 181 DO jj=1,jpj 182 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 183 END DO 184 185 DO jk = 1, jpkm1 186 DO jj = 1, jpj 187 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 188 & + uu_b(ji,jj,Krhs_a)-zub(ji,jj))*umask(ji,jj,jk) 189 END DO 190 END DO 191 END DO 192 193 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 194 ibdy1 = jpiglo-nbghostcells 195 ibdy2 = jpiglo-1 196 DO ji = mi0(ibdy1), mi1(ibdy2) 197 zvb(ji,:) = 0._wp 198 DO jk = 1, jpkm1 194 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 195 & + uu_b(ji,jj,Krhs_a)-zub(ji,jj))*umask(ji,jj,jk) 196 END DO 197 END DO 198 END DO 199 200 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 201 ibdy1 = jpiglo-nbghostcells 202 ibdy2 = jpiglo-1 203 DO ji = mi0(ibdy1), mi1(ibdy2) 204 zvb(ji,:) = 0._wp 205 DO jk = 1, jpkm1 206 DO jj = 1, jpj 207 zvb(ji,jj) = zvb(ji,jj) & 208 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 209 END DO 210 END DO 199 211 DO jj = 1, jpj 200 zvb(ji,jj) = zvb(ji,jj) & 212 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 213 END DO 214 DO jk = 1, jpkm1 215 DO jj = 1, jpj 216 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 217 & + vv_b(ji,jj,Krhs_a)-zvb(ji,jj)) * vmask(ji,jj,jk) 218 END DO 219 END DO 220 END DO 221 ENDIF 222 ENDIF 223 224 ! --- South --- ! 225 IF( lk_south ) THEN 226 jbdy1 = 2 227 jbdy2 = 1+nbghostcells 228 ! 229 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 230 DO jj = mj0(jbdy1), mj1(jbdy2) 231 vv_b(:,jj,Krhs_a) = 0._wp 232 DO jk = 1, jpkm1 233 DO ji = 1, jpi 234 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 235 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 236 END DO 237 END DO 238 DO ji=1,jpi 239 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 240 END DO 241 END DO 242 ENDIF 243 ! 244 DO jj = mj0(jbdy1), mj1(jbdy2) 245 zvb(:,jj) = 0._wp ! Correct transport 246 DO jk=1,jpkm1 247 DO ji=1,jpi 248 zvb(ji,jj) = zvb(ji,jj) & 201 249 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 202 250 END DO 203 251 END DO 204 DO j j = 1, jpj252 DO ji = 1, jpi 205 253 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 206 254 END DO 207 DO jk = 1, jpkm1 208 DO jj = 1, jpj 209 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 210 & + vv_b(ji,jj,Krhs_a)-zvb(ji,jj)) * vmask(ji,jj,jk) 211 END DO 212 END DO 213 END DO 214 ENDIF 215 216 ! --- South --- ! 217 jbdy1 = 2 218 jbdy2 = 1+nbghostcells 219 ! 220 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 221 DO jj = mj0(jbdy1), mj1(jbdy2) 222 vv_b(:,jj,Krhs_a) = 0._wp 255 223 256 DO jk = 1, jpkm1 224 257 DO ji = 1, jpi 225 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 226 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 227 END DO 228 END DO 229 DO ji=1,jpi 230 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 231 END DO 232 END DO 233 ENDIF 234 ! 235 DO jj = mj0(jbdy1), mj1(jbdy2) 236 zvb(:,jj) = 0._wp ! Correct transport 237 DO jk=1,jpkm1 238 DO ji=1,jpi 239 zvb(ji,jj) = zvb(ji,jj) & 240 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 241 END DO 242 END DO 243 DO ji = 1, jpi 244 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 245 END DO 246 247 DO jk = 1, jpkm1 258 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 259 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 260 END DO 261 END DO 262 END DO 263 264 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 265 DO jj = mj0(jbdy1), mj1(jbdy2) 266 zub(:,jj) = 0._wp 267 DO jk = 1, jpkm1 268 DO ji = 1, jpi 269 zub(ji,jj) = zub(ji,jj) & 270 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 271 END DO 272 END DO 273 DO ji = 1, jpi 274 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 275 END DO 276 277 DO jk = 1, jpkm1 278 DO ji = 1, jpi 279 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 280 & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 281 END DO 282 END DO 283 END DO 284 ENDIF 285 ENDIF 286 287 ! --- North --- ! 288 IF( lk_north ) THEN 289 jbdy1 = jpjglo-1-nbghostcells 290 jbdy2 = jpjglo-2 291 ! 292 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 293 DO jj = mj0(jbdy1), mj1(jbdy2) 294 vv_b(:,jj,Krhs_a) = 0._wp 295 DO jk = 1, jpkm1 296 DO ji = 1, jpi 297 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 298 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 299 END DO 300 END DO 301 DO ji=1,jpi 302 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 303 END DO 304 END DO 305 ENDIF 306 ! 307 DO jj = mj0(jbdy1), mj1(jbdy2) 308 zvb(:,jj) = 0._wp ! Correct transport 309 DO jk=1,jpkm1 310 DO ji=1,jpi 311 zvb(ji,jj) = zvb(ji,jj) & 312 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 313 END DO 314 END DO 248 315 DO ji = 1, jpi 249 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 250 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 251 END DO 252 END DO 253 END DO 254 255 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 256 DO jj = mj0(jbdy1), mj1(jbdy2) 257 zub(:,jj) = 0._wp 316 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 317 END DO 318 258 319 DO jk = 1, jpkm1 259 320 DO ji = 1, jpi 260 zub(ji,jj) = zub(ji,jj) & 261 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 262 END DO 263 END DO 264 DO ji = 1, jpi 265 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 266 END DO 321 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 322 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 323 END DO 324 END DO 325 END DO 267 326 268 DO jk = 1, jpkm1 327 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 328 jbdy1 = jpjglo-nbghostcells 329 jbdy2 = jpjglo-1 330 DO jj = mj0(jbdy1), mj1(jbdy2) 331 zub(:,jj) = 0._wp 332 DO jk = 1, jpkm1 333 DO ji = 1, jpi 334 zub(ji,jj) = zub(ji,jj) & 335 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 336 END DO 337 END DO 269 338 DO ji = 1, jpi 270 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 271 & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 272 END DO 273 END DO 274 END DO 275 ENDIF 276 277 ! --- North --- ! 278 jbdy1 = jpjglo-1-nbghostcells 279 jbdy2 = jpjglo-2 280 ! 281 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 282 DO jj = mj0(jbdy1), mj1(jbdy2) 283 vv_b(:,jj,Krhs_a) = 0._wp 284 DO jk = 1, jpkm1 285 DO ji = 1, jpi 286 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 287 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 288 END DO 289 END DO 290 DO ji=1,jpi 291 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 292 END DO 293 END DO 294 ENDIF 295 ! 296 DO jj = mj0(jbdy1), mj1(jbdy2) 297 zvb(:,jj) = 0._wp ! Correct transport 298 DO jk=1,jpkm1 299 DO ji=1,jpi 300 zvb(ji,jj) = zvb(ji,jj) & 301 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 302 END DO 303 END DO 304 DO ji = 1, jpi 305 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 306 END DO 307 308 DO jk = 1, jpkm1 309 DO ji = 1, jpi 310 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 311 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 312 END DO 313 END DO 314 END DO 315 316 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 317 jbdy1 = jpjglo-nbghostcells 318 jbdy2 = jpjglo-1 319 DO jj = mj0(jbdy1), mj1(jbdy2) 320 zub(:,jj) = 0._wp 321 DO jk = 1, jpkm1 322 DO ji = 1, jpi 323 zub(ji,jj) = zub(ji,jj) & 324 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 325 END DO 326 END DO 327 DO ji = 1, jpi 328 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 329 END DO 330 331 DO jk = 1, jpkm1 332 DO ji = 1, jpi 333 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 334 & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 335 END DO 336 END DO 337 END DO 339 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 340 END DO 341 342 DO jk = 1, jpkm1 343 DO ji = 1, jpi 344 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 345 & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 346 END DO 347 END DO 348 END DO 349 ENDIF 338 350 ENDIF 339 351 ! … … 354 366 ! 355 367 !--- West ---! 356 istart = 2 357 iend = nbghostcells+1 358 DO ji = mi0(istart), mi1(iend) 359 DO jj=1,jpj 360 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 361 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 362 END DO 363 END DO 368 IF( lk_west ) THEN 369 istart = 2 370 iend = nbghostcells+1 371 DO ji = mi0(istart), mi1(iend) 372 DO jj=1,jpj 373 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 374 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 375 END DO 376 END DO 377 ENDIF 364 378 ! 365 379 !--- East ---! 366 istart = jpiglo-nbghostcells 367 iend = jpiglo-1 368 DO ji = mi0(istart), mi1(iend) 369 DO jj=1,jpj 370 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 371 END DO 372 END DO 373 istart = jpiglo-nbghostcells-1 374 iend = jpiglo-2 375 DO ji = mi0(istart), mi1(iend) 376 DO jj=1,jpj 377 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 378 END DO 379 END DO 380 IF( lk_east ) THEN 381 istart = jpiglo-nbghostcells 382 iend = jpiglo-1 383 DO ji = mi0(istart), mi1(iend) 384 385 DO jj=1,jpj 386 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 387 END DO 388 END DO 389 istart = jpiglo-nbghostcells-1 390 iend = jpiglo-2 391 DO ji = mi0(istart), mi1(iend) 392 DO jj=1,jpj 393 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 394 END DO 395 END DO 396 ENDIF 380 397 ! 381 398 !--- South ---! 382 jstart = 2 383 jend = nbghostcells+1 384 DO jj = mj0(jstart), mj1(jend) 385 DO ji=1,jpi 386 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 387 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 388 END DO 389 END DO 399 IF( lk_south ) THEN 400 jstart = 2 401 jend = nbghostcells+1 402 DO jj = mj0(jstart), mj1(jend) 403 404 DO ji=1,jpi 405 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 406 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 407 END DO 408 END DO 409 ENDIF 390 410 ! 391 411 !--- North ---! 392 jstart = jpjglo-nbghostcells 393 jend = jpjglo-1 394 DO jj = mj0(jstart), mj1(jend) 395 DO ji=1,jpi 396 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 397 END DO 398 END DO 399 jstart = jpjglo-nbghostcells-1 400 jend = jpjglo-2 401 DO jj = mj0(jstart), mj1(jend) 402 DO ji=1,jpi 403 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 404 END DO 405 END DO 412 IF( lk_north ) THEN 413 jstart = jpjglo-nbghostcells 414 jend = jpjglo-1 415 DO jj = mj0(jstart), mj1(jend) 416 DO ji=1,jpi 417 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 418 END DO 419 END DO 420 jstart = jpjglo-nbghostcells-1 421 jend = jpjglo-2 422 DO jj = mj0(jstart), mj1(jend) 423 DO ji=1,jpi 424 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 425 END DO 426 END DO 427 ENDIF 406 428 ! 407 429 END SUBROUTINE Agrif_dyn_ts … … 421 443 ! 422 444 !--- West ---! 423 istart = 2 424 iend = nbghostcells+1 425 DO ji = mi0(istart), mi1(iend) 426 DO jj=1,jpj 427 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 428 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 429 END DO 430 END DO 445 IF( lk_west ) THEN 446 istart = 2 447 iend = nbghostcells+1 448 DO ji = mi0(istart), mi1(iend) 449 DO jj=1,jpj 450 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 451 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 452 END DO 453 END DO 454 ENDIF 431 455 ! 432 456 !--- East ---! 433 istart = jpiglo-nbghostcells 434 iend = jpiglo-1 435 DO ji = mi0(istart), mi1(iend) 436 DO jj=1,jpj 437 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 438 END DO 439 END DO 440 istart = jpiglo-nbghostcells-1 441 iend = jpiglo-2 442 DO ji = mi0(istart), mi1(iend) 443 DO jj=1,jpj 444 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 445 END DO 446 END DO 457 IF( lk_east ) THEN 458 istart = jpiglo-nbghostcells 459 iend = jpiglo-1 460 DO ji = mi0(istart), mi1(iend) 461 DO jj=1,jpj 462 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 463 END DO 464 END DO 465 istart = jpiglo-nbghostcells-1 466 iend = jpiglo-2 467 DO ji = mi0(istart), mi1(iend) 468 DO jj=1,jpj 469 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 470 END DO 471 END DO 472 ENDIF 447 473 ! 448 474 !--- South ---! 449 jstart = 2 450 jend = nbghostcells+1 451 DO jj = mj0(jstart), mj1(jend) 452 DO ji=1,jpi 453 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 454 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 455 END DO 456 END DO 475 IF( lk_south ) THEN 476 jstart = 2 477 jend = nbghostcells+1 478 DO jj = mj0(jstart), mj1(jend) 479 DO ji=1,jpi 480 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 481 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 482 END DO 483 END DO 484 ENDIF 457 485 ! 458 486 !--- North ---! 459 jstart = jpjglo-nbghostcells 460 jend = jpjglo-1 461 DO jj = mj0(jstart), mj1(jend) 462 DO ji=1,jpi 463 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 464 END DO 465 END DO 466 jstart = jpjglo-nbghostcells-1 467 jend = jpjglo-2 468 DO jj = mj0(jstart), mj1(jend) 469 DO ji=1,jpi 470 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 471 END DO 472 END DO 487 IF( lk_north ) THEN 488 jstart = jpjglo-nbghostcells 489 jend = jpjglo-1 490 DO jj = mj0(jstart), mj1(jend) 491 DO ji=1,jpi 492 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 493 END DO 494 END DO 495 jstart = jpjglo-nbghostcells-1 496 jend = jpjglo-2 497 DO jj = mj0(jstart), mj1(jend) 498 DO ji=1,jpi 499 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 500 END DO 501 END DO 502 ENDIF 473 503 ! 474 504 END SUBROUTINE Agrif_dyn_ts_flux … … 494 524 Agrif_SpecialValue = 0._wp 495 525 Agrif_UseSpecialValue = ln_spc_dyn 526 527 use_sign_north = .TRUE. 528 sign_north = -1. 529 496 530 ! 497 531 ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners) … … 518 552 ENDIF 519 553 Agrif_UseSpecialValue = .FALSE. 554 use_sign_north = .FALSE. 520 555 ! 521 556 END SUBROUTINE Agrif_dta_ts … … 542 577 ! 543 578 ! --- West --- ! 544 istart = 2 545 iend = 1 + nbghostcells 546 DO ji = mi0(istart), mi1(iend) 547 DO jj = 1, jpj 548 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 579 IF(lk_west) THEN 580 istart = 2 581 iend = 1 + nbghostcells 582 DO ji = mi0(istart), mi1(iend) 583 DO jj = 1, jpj 584 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 585 ENDDO 549 586 ENDDO 550 END DO587 ENDIF 551 588 ! 552 589 ! --- East --- ! 553 istart = jpiglo - nbghostcells 554 iend = jpiglo - 1 555 DO ji = mi0(istart), mi1(iend) 556 DO jj = 1, jpj 557 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 590 IF(lk_east) THEN 591 istart = jpiglo - nbghostcells 592 iend = jpiglo - 1 593 DO ji = mi0(istart), mi1(iend) 594 DO jj = 1, jpj 595 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 596 ENDDO 558 597 ENDDO 559 END DO598 ENDIF 560 599 ! 561 600 ! --- South --- ! 562 jstart = 2 563 jend = 1 + nbghostcells 564 DO jj = mj0(jstart), mj1(jend) 565 DO ji = 1, jpi 566 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 601 IF(lk_south) THEN 602 jstart = 2 603 jend = 1 + nbghostcells 604 DO jj = mj0(jstart), mj1(jend) 605 DO ji = 1, jpi 606 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 607 ENDDO 567 608 ENDDO 568 END DO609 ENDIF 569 610 ! 570 611 ! --- North --- ! 571 jstart = jpjglo - nbghostcells 572 jend = jpjglo - 1 573 DO jj = mj0(jstart), mj1(jend) 574 DO ji = 1, jpi 575 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 612 IF(lk_north) THEN 613 jstart = jpjglo - nbghostcells 614 jend = jpjglo - 1 615 DO jj = mj0(jstart), mj1(jend) 616 DO ji = 1, jpi 617 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 618 ENDDO 576 619 ENDDO 577 END DO620 ENDIF 578 621 ! 579 622 END SUBROUTINE Agrif_ssh … … 593 636 ! 594 637 ! --- West --- ! 595 istart = 2 596 iend = 1+nbghostcells 597 DO ji = mi0(istart), mi1(iend) 598 DO jj = 1, jpj 599 ssha_e(ji,jj) = hbdy(ji,jj) 638 IF(lk_west) THEN 639 istart = 2 640 iend = 1+nbghostcells 641 DO ji = mi0(istart), mi1(iend) 642 DO jj = 1, jpj 643 ssha_e(ji,jj) = hbdy(ji,jj) 644 ENDDO 600 645 ENDDO 601 END DO646 ENDIF 602 647 ! 603 648 ! --- East --- ! 604 istart = jpiglo - nbghostcells 605 iend = jpiglo - 1 606 DO ji = mi0(istart), mi1(iend) 607 DO jj = 1, jpj 608 ssha_e(ji,jj) = hbdy(ji,jj) 649 IF(lk_east) THEN 650 istart = jpiglo - nbghostcells 651 iend = jpiglo - 1 652 DO ji = mi0(istart), mi1(iend) 653 DO jj = 1, jpj 654 ssha_e(ji,jj) = hbdy(ji,jj) 655 ENDDO 609 656 ENDDO 610 END DO657 ENDIF 611 658 ! 612 659 ! --- South --- ! 613 jstart = 2 614 jend = 1+nbghostcells 615 DO jj = mj0(jstart), mj1(jend) 616 DO ji = 1, jpi 617 ssha_e(ji,jj) = hbdy(ji,jj) 660 IF(lk_south) THEN 661 jstart = 2 662 jend = 1+nbghostcells 663 DO jj = mj0(jstart), mj1(jend) 664 DO ji = 1, jpi 665 ssha_e(ji,jj) = hbdy(ji,jj) 666 ENDDO 618 667 ENDDO 619 END DO668 ENDIF 620 669 ! 621 670 ! --- North --- ! 622 jstart = jpjglo - nbghostcells 623 jend = jpjglo - 1 624 DO jj = mj0(jstart), mj1(jend) 625 DO ji = 1, jpi 626 ssha_e(ji,jj) = hbdy(ji,jj) 671 IF(lk_north) THEN 672 jstart = jpjglo - nbghostcells 673 jend = jpjglo - 1 674 DO jj = mj0(jstart), mj1(jend) 675 DO ji = 1, jpi 676 ssha_e(ji,jj) = hbdy(ji,jj) 677 ENDDO 627 678 ENDDO 628 END DO679 ENDIF 629 680 ! 630 681 END SUBROUTINE Agrif_ssh_ts … … 662 713 INTEGER :: ji, jj, jk, jn ! dummy loop indices 663 714 INTEGER :: N_in, N_out 715 INTEGER :: item 664 716 ! vertical interpolation: 665 717 REAL(wp) :: zhtot 666 718 REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin 667 REAL(wp), DIMENSION(k1:k2) :: h_in 668 REAL(wp), DIMENSION(1:jpk) :: h_out 669 !!---------------------------------------------------------------------- 670 671 IF( before ) THEN 719 REAL(wp), DIMENSION(k1:k2) :: h_in, z_in 720 REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 721 !!---------------------------------------------------------------------- 722 723 IF( before ) THEN 724 725 item = Kmm_a 726 IF( l_ini_child ) Kmm_a = Kbb_a 727 672 728 DO jn = 1,jpts 673 729 DO jk=k1,k2 … … 678 734 END DO 679 735 END DO 680 END DO 681 682 # if defined key_vertical 683 ! Interpolate thicknesses 684 ! Warning: these are masked, hence extrapolated prior interpolation. 685 DO jk=k1,k2 686 DO jj=j1,j2 687 DO ji=i1,i2 688 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 689 END DO 690 END DO 691 END DO 692 693 ! Extrapolate thicknesses in partial bottom cells: 694 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 695 IF (ln_zps) THEN 696 DO jj=j1,j2 697 DO ji=i1,i2 698 jk = mbkt(ji,jj) 699 ptab(ji,jj,jk,jpts+1) = 0._wp 700 END DO 701 END DO 702 END IF 703 704 ! Save ssh at last level: 705 IF (.NOT.ln_linssh) THEN 706 ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 707 ELSE 708 ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 709 END IF 710 # endif 736 END DO 737 738 IF( l_vremap .OR. l_ini_child) THEN 739 ! Interpolate thicknesses 740 ! Warning: these are masked, hence extrapolated prior interpolation. 741 DO jk=k1,k2 742 DO jj=j1,j2 743 DO ji=i1,i2 744 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 745 746 END DO 747 END DO 748 END DO 749 750 ! Extrapolate thicknesses in partial bottom cells: 751 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 752 IF (ln_zps) THEN 753 DO jj=j1,j2 754 DO ji=i1,i2 755 jk = mbkt(ji,jj) 756 ptab(ji,jj,jk,jpts+1) = 0._wp 757 END DO 758 END DO 759 END IF 760 761 ! Save ssh at last level: 762 IF (.NOT.ln_linssh) THEN 763 ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 764 ELSE 765 ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 766 END IF 767 ENDIF 768 Kmm_a = item 769 711 770 ELSE 712 713 # if defined key_vertical 714 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 715 716 DO jj=j1,j2 717 DO ji=i1,i2 718 ts(ji,jj,:,:,Krhs_a) = 0._wp 719 N_in = mbkt_parent(ji,jj) 720 zhtot = 0._wp 721 DO jk=1,N_in !k2 = jpk of parent grid 722 IF (jk==N_in) THEN 723 h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 724 ELSE 725 h_in(jk) = ptab(ji,jj,jk,n2) 771 item = Krhs_a 772 IF( l_ini_child ) Krhs_a = Kbb_a 773 774 IF( l_vremap .OR. l_ini_child ) THEN 775 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 776 777 DO jj=j1,j2 778 DO ji=i1,i2 779 ts(ji,jj,:,:,Krhs_a) = 0. 780 ! IF( l_ini_child) ts(ji,jj,:,:,Krhs_a) = ptab(ji,jj,:,1:jpts) 781 N_in = mbkt_parent(ji,jj) 782 zhtot = 0._wp 783 DO jk=1,N_in !k2 = jpk of parent grid 784 IF (jk==N_in) THEN 785 h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 786 ELSE 787 h_in(jk) = ptab(ji,jj,jk,n2) 788 ENDIF 789 zhtot = zhtot + h_in(jk) 790 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 791 END DO 792 z_in(1) = 0.5_wp * h_in(1) - zhtot + ht0_parent(ji,jj) 793 DO jk=2,N_in 794 z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 795 ENDDO 796 797 N_out = 0 798 DO jk=1,jpk ! jpk of child grid 799 IF (tmask(ji,jj,jk) == 0._wp) EXIT 800 N_out = N_out + 1 801 h_out(jk) = e3t(ji,jj,jk,Krhs_a) 802 ENDDO 803 804 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + ht_0(ji,jj) 805 DO jk=2,N_out 806 z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 807 ENDDO 808 809 IF (N_in*N_out > 0) THEN 810 IF( l_ini_child ) THEN 811 CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a), & 812 & z_out(1:N_out),N_in,N_out,jpts) 813 ELSE 814 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a), & 815 & h_out(1:N_out),N_in,N_out,jpts) 816 ENDIF 726 817 ENDIF 727 zhtot = zhtot + h_in(jk)728 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)729 END DO730 N_out = 0731 DO jk=1,jpk ! jpk of child grid732 IF (tmask(ji,jj,jk) == 0._wp) EXIT733 N_out = N_out + 1734 h_out(jk) = e3t(ji,jj,jk,Krhs_a)735 818 ENDDO 736 IF (N_in*N_out > 0) THEN737 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),h_out(1:N_out),N_in,N_out,jpts)738 ENDIF739 819 ENDDO 740 ENDDO 741 # else 742 ! 743 DO jn=1, jpts 744 ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 745 END DO 746 # endif 820 Krhs_a = item 821 822 ELSE 823 824 DO jn=1, jpts 825 ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 826 END DO 827 ENDIF 747 828 748 829 ENDIF … … 780 861 REAL(wp) :: zrhoy, zhtot 781 862 ! vertical interpolation: 782 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 783 REAL(wp), DIMENSION(1:jpk) :: h_out 784 INTEGER :: N_in, N_out 863 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in 864 REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 865 INTEGER :: N_in, N_out,item 785 866 REAL(wp) :: h_diff 786 867 !!--------------------------------------------- 787 868 ! 788 869 IF (before) THEN 870 871 item = Kmm_a 872 IF( l_ini_child ) Kmm_a = Kbb_a 873 789 874 DO jk=1,jpk 790 875 DO jj=j1,j2 791 876 DO ji=i1,i2 792 877 ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 793 # if defined key_vertical 794 ! Interpolate thicknesses (masked for subsequent extrapolation) 795 ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 796 # endif 797 END DO 798 END DO 799 END DO 800 # if defined key_vertical 878 IF( l_vremap .OR. l_ini_child) THEN 879 ! Interpolate thicknesses (masked for subsequent extrapolation) 880 ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 881 ENDIF 882 END DO 883 END DO 884 END DO 885 886 IF( l_vremap .OR. l_ini_child) THEN 801 887 ! Extrapolate thicknesses in partial bottom cells: 802 888 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 803 IF (ln_zps) THEN 804 DO jj=j1,j2 805 DO ji=i1,i2 806 jk = mbku(ji,jj) 807 ptab(ji,jj,jk,2) = 0._wp 808 END DO 809 END DO 810 END IF 811 ! Save ssh at last level: 812 ptab(i1:i2,j1:j2,k2,2) = 0._wp 813 IF (.NOT.ln_linssh) THEN 814 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 815 DO jk=1,jpk 816 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk) 817 END DO 818 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2) 819 END IF 820 # endif 889 IF (ln_zps) THEN 890 DO jj=j1,j2 891 DO ji=i1,i2 892 jk = mbku(ji,jj) 893 ptab(ji,jj,jk,2) = 0._wp 894 END DO 895 END DO 896 END IF 897 898 ! Save ssh at last level: 899 ptab(i1:i2,j1:j2,k2,2) = 0._wp 900 IF (.NOT.ln_linssh) THEN 901 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 902 DO jk=1,jpk 903 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk) 904 END DO 905 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2) 906 END IF 907 ENDIF 908 909 Kmm_a = item 821 910 ! 822 911 ELSE 823 912 zrhoy = Agrif_rhoy() 824 # if defined key_vertical 913 914 IF( l_vremap .OR. l_ini_child) THEN 825 915 ! VERTICAL REFINEMENT BEGIN 826 916 827 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 828 829 DO ji=i1,i2 830 DO jj=j1,j2 831 uu(ji,jj,:,Krhs_a) = 0._wp 832 N_in = mbku_parent(ji,jj) 833 zhtot = 0._wp 834 DO jk=1,N_in 835 IF (jk==N_in) THEN 836 h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 837 ELSE 838 h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 839 ENDIF 840 zhtot = zhtot + h_in(jk) 841 tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk)) 842 ENDDO 843 844 N_out = 0 845 DO jk=1,jpk 846 if (umask(ji,jj,jk) == 0) EXIT 847 N_out = N_out + 1 848 h_out(N_out) = e3u(ji,jj,jk,Krhs_a) 849 ENDDO 850 IF (N_in*N_out > 0) THEN 851 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 852 ENDIF 917 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 918 919 DO ji=i1,i2 920 DO jj=j1,j2 921 uu(ji,jj,:,Krhs_a) = 0._wp 922 N_in = mbku_parent(ji,jj) 923 zhtot = 0._wp 924 DO jk=1,N_in 925 IF (jk==N_in) THEN 926 h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 927 ELSE 928 h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 929 ENDIF 930 zhtot = zhtot + h_in(jk) 931 IF( h_in(jk) .GT. 0. ) THEN 932 tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk)) 933 ELSE 934 tabin(jk) = 0. 935 ENDIF 936 ENDDO 937 z_in(1) = 0.5_wp * h_in(1) - zhtot + hu0_parent(ji,jj) 938 DO jk=2,N_in 939 z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 940 ENDDO 941 942 N_out = 0 943 DO jk=1,jpk 944 IF (umask(ji,jj,jk) == 0) EXIT 945 N_out = N_out + 1 946 h_out(N_out) = e3u(ji,jj,jk,Krhs_a) 947 ENDDO 948 949 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hu_0(ji,jj) 950 DO jk=2,N_out 951 z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 952 ENDDO 953 954 IF (N_in*N_out > 0) THEN 955 IF( l_ini_child ) THEN 956 CALL remap_linear (tabin(1:N_in),z_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1) 957 ELSE 958 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 959 ENDIF 960 ENDIF 961 ENDDO 853 962 ENDDO 854 ENDDO 855 856 # else 857 DO jk = 1, jpkm1 858 DO jj=j1,j2 859 uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) ) 860 END DO 861 END DO 862 # endif 963 ELSE 964 DO jk = 1, jpkm1 965 DO jj=j1,j2 966 uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) ) 967 END DO 968 END DO 969 ENDIF 863 970 864 971 ENDIF … … 878 985 REAL(wp) :: zrhox 879 986 ! vertical interpolation: 880 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 881 REAL(wp), DIMENSION(1:jpk) :: h_out 882 INTEGER :: N_in, N_out 987 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in 988 REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 989 INTEGER :: N_in, N_out, item 883 990 REAL(wp) :: h_diff, zhtot 884 991 !!--------------------------------------------- 885 992 ! 886 IF (before) THEN 993 IF (before) THEN 994 995 item = Kmm_a 996 IF( l_ini_child ) Kmm_a = Kbb_a 997 887 998 DO jk=k1,k2 888 999 DO jj=j1,j2 889 1000 DO ji=i1,i2 890 1001 ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk)) 891 # if defined key_vertical 892 ! Interpolate thicknesses (masked for subsequent extrapolation) 893 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 894 # endif 895 END DO 896 END DO 897 END DO 898 # if defined key_vertical 1002 IF( l_vremap .OR. l_ini_child) THEN 1003 ! Interpolate thicknesses (masked for subsequent extrapolation) 1004 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 1005 ENDIF 1006 END DO 1007 END DO 1008 END DO 1009 1010 IF( l_vremap .OR. l_ini_child) THEN 899 1011 ! Extrapolate thicknesses in partial bottom cells: 900 1012 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 901 IF (ln_zps) THEN 1013 IF (ln_zps) THEN 1014 DO jj=j1,j2 1015 DO ji=i1,i2 1016 jk = mbkv(ji,jj) 1017 ptab(ji,jj,jk,2) = 0._wp 1018 END DO 1019 END DO 1020 END IF 1021 ! Save ssh at last level: 1022 ptab(i1:i2,j1:j2,k2,2) = 0._wp 1023 IF (.NOT.ln_linssh) THEN 1024 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 1025 DO jk=1,jpk 1026 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 1027 END DO 1028 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2) 1029 END IF 1030 ENDIF 1031 item = Kmm_a 1032 1033 ELSE 1034 zrhox = Agrif_rhox() 1035 1036 IF( l_vremap .OR. l_ini_child ) THEN 1037 1038 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1039 902 1040 DO jj=j1,j2 903 1041 DO ji=i1,i2 904 jk = mbkv(ji,jj) 905 ptab(ji,jj,jk,2) = 0._wp 906 END DO 907 END DO 908 END IF 909 ! Save ssh at last level: 910 ptab(i1:i2,j1:j2,k2,2) = 0._wp 911 IF (.NOT.ln_linssh) THEN 912 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 913 DO jk=1,jpk 914 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 915 END DO 916 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2) 917 END IF 918 # endif 919 ELSE 920 zrhox = Agrif_rhox() 921 # if defined key_vertical 922 923 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 924 925 DO jj=j1,j2 926 DO ji=i1,i2 927 vv(ji,jj,:,Krhs_a) = 0._wp 928 N_in = mbkv_parent(ji,jj) 929 zhtot = 0._wp 930 DO jk=1,N_in 931 IF (jk==N_in) THEN 932 h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 933 ELSE 934 h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 1042 vv(ji,jj,:,Krhs_a) = 0._wp 1043 N_in = mbkv_parent(ji,jj) 1044 zhtot = 0._wp 1045 DO jk=1,N_in 1046 IF (jk==N_in) THEN 1047 h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 1048 ELSE 1049 h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 1050 ENDIF 1051 zhtot = zhtot + h_in(jk) 1052 IF( h_in(jk) .GT. 0. ) THEN 1053 tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk)) 1054 ELSE 1055 tabin(jk) = 0. 1056 ENDIF 1057 ENDDO 1058 1059 z_in(1) = 0.5_wp * h_in(1) - zhtot + hv0_parent(ji,jj) 1060 DO jk=2,N_in 1061 z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 1062 ENDDO 1063 1064 N_out = 0 1065 DO jk=1,jpk 1066 IF (vmask(ji,jj,jk) == 0) EXIT 1067 N_out = N_out + 1 1068 h_out(N_out) = e3v(ji,jj,jk,Krhs_a) 1069 ENDDO 1070 1071 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hv_0(ji,jj) 1072 DO jk=2,N_out 1073 z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 1074 ENDDO 1075 1076 IF (N_in*N_out > 0) THEN 1077 IF( l_ini_child ) THEN 1078 CALL remap_linear (tabin(1:N_in),z_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1) 1079 ELSE 1080 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 1081 ENDIF 935 1082 ENDIF 936 zhtot = zhtot + h_in(jk) 937 tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk)) 938 ENDDO 939 940 N_out = 0 941 DO jk=1,jpk 942 if (vmask(ji,jj,jk) == 0) EXIT 943 N_out = N_out + 1 944 h_out(N_out) = e3v(ji,jj,jk,Krhs_a) 945 END DO 946 IF (N_in*N_out > 0) THEN 947 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 948 ENDIF 949 END DO 950 END DO 951 # else 952 DO jk = 1, jpkm1 953 vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) ) 954 END DO 955 # endif 1083 END DO 1084 END DO 1085 ELSE 1086 DO jk = 1, jpkm1 1087 vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) ) 1088 END DO 1089 ENDIF 956 1090 ENDIF 957 1091 ! … … 1163 1297 END SUBROUTINE interpe3t 1164 1298 1165 1166 1299 SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 1167 1300 !!---------------------------------------------------------------------- … … 1185 1318 END DO 1186 1319 END DO 1187 END DO 1188 1189 # if defined key_vertical 1190 ! Interpolate thicknesses 1191 ! Warning: these are masked, hence extrapolated prior interpolation. 1192 DO jk=k1,k2 1193 DO jj=j1,j2 1194 DO ji=i1,i2 1195 ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 1196 END DO 1197 END DO 1198 END DO 1199 1200 ! Extrapolate thicknesses in partial bottom cells: 1201 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 1202 IF (ln_zps) THEN 1203 DO jj=j1,j2 1204 DO ji=i1,i2 1205 jk = mbkt(ji,jj) 1206 ptab(ji,jj,jk,2) = 0._wp 1207 END DO 1208 END DO 1209 END IF 1210 1211 ! Save ssh at last level: 1212 IF (.NOT.ln_linssh) THEN 1213 ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 1214 ELSE 1215 ptab(i1:i2,j1:j2,k2,2) = 0._wp 1216 END IF 1217 # endif 1320 END DO 1321 1322 IF( l_vremap ) THEN 1323 ! Interpolate thicknesses 1324 ! Warning: these are masked, hence extrapolated prior interpolation. 1325 DO jk=k1,k2 1326 DO jj=j1,j2 1327 DO ji=i1,i2 1328 ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 1329 END DO 1330 END DO 1331 END DO 1332 1333 ! Extrapolate thicknesses in partial bottom cells: 1334 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 1335 IF (ln_zps) THEN 1336 DO jj=j1,j2 1337 DO ji=i1,i2 1338 jk = mbkt(ji,jj) 1339 ptab(ji,jj,jk,2) = 0._wp 1340 END DO 1341 END DO 1342 END IF 1343 1344 ! Save ssh at last level: 1345 IF (.NOT.ln_linssh) THEN 1346 ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 1347 ELSE 1348 ptab(i1:i2,j1:j2,k2,2) = 0._wp 1349 END IF 1350 ENDIF 1351 1218 1352 ELSE 1219 #ifdef key_vertical 1220 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1221 avm_k(i1:i2,j1:j2,k1:k2) = 0._wp 1222 1223 DO jj = j1, j2 1224 DO ji =i1, i2 1225 N_in = mbkt_parent(ji,jj) 1226 IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 1227 z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2) 1228 DO jk = N_in, 1, -1 ! Parent vertical grid 1229 z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2) 1230 tabin(jk) = ptab(ji,jj,jk,1) 1231 END DO 1232 N_out = mbkt(ji,jj) 1233 DO jk = 1, N_out ! Child vertical grid 1234 z_out(jk) = gdepw(ji,jj,jk,Kmm_a) 1353 1354 IF( l_vremap ) THEN 1355 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1356 avm_k(i1:i2,j1:j2,k1:k2) = 0._wp 1357 1358 DO jj = j1, j2 1359 DO ji =i1, i2 1360 N_in = mbkt_parent(ji,jj) 1361 IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 1362 z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2) 1363 DO jk = N_in, 1, -1 ! Parent vertical grid 1364 z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2) 1365 tabin(jk) = ptab(ji,jj,jk,1) 1366 END DO 1367 N_out = mbkt(ji,jj) 1368 DO jk = 1, N_out ! Child vertical grid 1369 z_out(jk) = gdepw(ji,jj,jk,Kmm_a) 1370 ENDDO 1371 IF (N_in*N_out > 0) THEN 1372 CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1) 1373 ENDIF 1235 1374 ENDDO 1236 IF (N_in*N_out > 0) THEN1237 CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1)1238 ENDIF1239 1375 ENDDO 1240 ENDDO 1241 #else 1242 avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1) 1243 #endif 1376 ELSE 1377 avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1) 1378 ENDIF 1244 1379 ENDIF 1245 1380 ! 1246 1381 END SUBROUTINE interpavm 1247 1382 1248 # if defined key_vertical1249 1383 SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before ) 1250 1384 !!---------------------------------------------------------------------- … … 1282 1416 ! 1283 1417 END SUBROUTINE interpht0 1284 #endif 1285 1418 1419 SUBROUTINE agrif_initts(tabres,i1,i2,j1,j2,k1,k2,m1,m2,before) 1420 INTEGER :: i1, i2, j1, j2, k1, k2, m1, m2 1421 REAL(wp):: tabres(i1:i2,j1:j2,k1:k2,m1:m2) 1422 LOGICAL :: before 1423 1424 INTEGER :: jm 1425 1426 IF (before) THEN 1427 DO jm=1,jpts 1428 tabres(i1:i2,j1:j2,k1:k2,jm) = ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a) 1429 END DO 1430 ELSE 1431 DO jm=1,jpts 1432 ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)=tabres(i1:i2,j1:j2,k1:k2,jm) 1433 END DO 1434 ENDIF 1435 END SUBROUTINE agrif_initts 1436 1437 SUBROUTINE agrif_initssh( ptab, i1, i2, j1, j2, before ) 1438 !!---------------------------------------------------------------------- 1439 !! *** ROUTINE interpsshn *** 1440 !!---------------------------------------------------------------------- 1441 INTEGER , INTENT(in ) :: i1, i2, j1, j2 1442 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1443 LOGICAL , INTENT(in ) :: before 1444 ! 1445 !!---------------------------------------------------------------------- 1446 ! 1447 IF( before) THEN 1448 ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kbb_a) 1449 ELSE 1450 ssh(i1:i2,j1:j2,Kbb_a) = ptab(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 1451 ENDIF 1452 ! 1453 END SUBROUTINE agrif_initssh 1454 1286 1455 #else 1287 1456 !!---------------------------------------------------------------------- -
NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/NST/agrif_oce_sponge.F90
r12546 r13220 80 80 Agrif_SpecialValue=0. 81 81 Agrif_UseSpecialValue = ln_spc_dyn 82 use_sign_north = .TRUE. 83 sign_north = -1. 82 84 ! 83 85 tabspongedone_u = .FALSE. … … 90 92 ! 91 93 Agrif_UseSpecialValue = .FALSE. 94 use_sign_north = .FALSE. 92 95 #endif 93 96 ! … … 127 130 128 131 ! --- West --- ! 129 ztabramp(:,:) = 0._wp 130 ind1 = 1+nbghostcells 131 DO ji = mi0(ind1), mi1(ind1) 132 ztabramp(ji,:) = ssumask(ji,:) 133 END DO 134 ! 135 zmskwest(:) = 0._wp 136 zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 132 IF( lk_west) THEN 133 ztabramp(:,:) = 0._wp 134 ind1 = 1+nbghostcells 135 DO ji = mi0(ind1), mi1(ind1) 136 ztabramp(ji,:) = ssumask(ji,:) 137 END DO 138 ! 139 zmskwest(:) = 0._wp 140 zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 141 ENDIF 137 142 138 143 ! --- East --- ! 139 ztabramp(:,:) = 0._wp 140 ind1 = jpiglo - nbghostcells - 1 141 DO ji = mi0(ind1), mi1(ind1) 142 ztabramp(ji,:) = ssumask(ji,:) 143 END DO 144 ! 145 zmskeast(:) = 0._wp 146 zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 144 IF( lk_east ) THEN 145 ztabramp(:,:) = 0._wp 146 ind1 = jpiglo - nbghostcells - 1 147 DO ji = mi0(ind1), mi1(ind1) 148 ztabramp(ji,:) = ssumask(ji,:) 149 END DO 150 ! 151 zmskeast(:) = 0._wp 152 zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 153 ENDIF 147 154 148 155 ! --- South --- ! 149 ztabramp(:,:) = 0._wp 150 ind1 = 1+nbghostcells 151 DO jj = mj0(ind1), mj1(ind1) 152 ztabramp(:,jj) = ssvmask(:,jj) 153 END DO 154 ! 155 zmsksouth(:) = 0._wp 156 zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 156 IF( lk_south ) THEN 157 ztabramp(:,:) = 0._wp 158 ind1 = 1+nbghostcells 159 DO jj = mj0(ind1), mj1(ind1) 160 ztabramp(:,jj) = ssvmask(:,jj) 161 END DO 162 ! 163 zmsksouth(:) = 0._wp 164 zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 165 ENDIF 157 166 158 167 ! --- North --- ! 159 ztabramp(:,:) = 0._wp 160 ind1 = jpjglo - nbghostcells - 1 161 DO jj = mj0(ind1), mj1(ind1) 162 ztabramp(:,jj) = ssvmask(:,jj) 163 END DO 164 ! 165 zmsknorth(:) = 0._wp 166 zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 168 IF( lk_north) THEN 169 ztabramp(:,:) = 0._wp 170 ind1 = jpjglo - nbghostcells - 1 171 DO jj = mj0(ind1), mj1(ind1) 172 ztabramp(:,jj) = ssvmask(:,jj) 173 END DO 174 ! 175 zmsknorth(:) = 0._wp 176 zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 177 ENDIF 178 167 179 ! JC: SPONGE MASKING TO BE SORTED OUT: 168 180 zmskwest(:) = 1._wp … … 192 204 193 205 ! --- West --- ! 194 ind1 = 1+nbghostcells 195 ind2 = 1+nbghostcells + ispongearea 196 DO ji = mi0(ind1), mi1(ind2) 197 DO jj = 1, jpj 198 ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj) 199 END DO 200 END DO 201 202 ! ghost cells: 203 ind1 = 1 204 ind2 = nbghostcells + 1 205 DO ji = mi0(ind1), mi1(ind2) 206 DO jj = 1, jpj 207 ztabramp(ji,jj) = zmskwest(jj) 208 END DO 209 END DO 206 IF(lk_west) THEN 207 ind1 = 1+nbghostcells 208 ind2 = 1+nbghostcells + ispongearea 209 DO ji = mi0(ind1), mi1(ind2) 210 DO jj = 1, jpj 211 ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj) 212 END DO 213 END DO 214 215 ! ghost cells: 216 ind1 = 1 217 ind2 = nbghostcells + 1 218 DO ji = mi0(ind1), mi1(ind2) 219 DO jj = 1, jpj 220 ztabramp(ji,jj) = zmskwest(jj) 221 END DO 222 END DO 223 ENDIF 210 224 211 225 ! --- East --- ! 212 ind1 = jpiglo - nbghostcells - ispongearea 213 ind2 = jpiglo - nbghostcells 214 DO ji = mi0(ind1), mi1(ind2) 215 DO jj = 1, jpj 216 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj) 217 ENDDO 218 END DO 219 220 ! ghost cells: 221 ind1 = jpiglo - nbghostcells 222 ind2 = jpiglo 223 DO ji = mi0(ind1), mi1(ind2) 224 DO jj = 1, jpj 225 ztabramp(ji,jj) = zmskeast(jj) 226 ENDDO 227 END DO 226 IF(lk_east) THEN 227 ind1 = jpiglo - nbghostcells - ispongearea 228 ind2 = jpiglo - nbghostcells 229 DO ji = mi0(ind1), mi1(ind2) 230 231 DO jj = 1, jpj 232 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj) 233 ENDDO 234 END DO 235 236 ! ghost cells: 237 ind1 = jpiglo - nbghostcells 238 ind2 = jpiglo 239 DO ji = mi0(ind1), mi1(ind2) 240 241 DO jj = 1, jpj 242 ztabramp(ji,jj) = zmskeast(jj) 243 ENDDO 244 END DO 245 ENDIF 228 246 229 247 ! --- South --- ! 230 ind1 = 1+nbghostcells 231 ind2 = 1+nbghostcells + jspongearea 232 DO jj = mj0(ind1), mj1(ind2) 233 DO ji = 1, jpi 234 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji) 235 END DO 236 END DO 237 238 ! ghost cells: 239 ind1 = 1 240 ind2 = nbghostcells + 1 241 DO jj = mj0(ind1), mj1(ind2) 242 DO ji = 1, jpi 243 ztabramp(ji,jj) = zmsksouth(ji) 244 END DO 245 END DO 248 IF( lk_south ) THEN 249 ind1 = 1+nbghostcells 250 ind2 = 1+nbghostcells + jspongearea 251 DO jj = mj0(ind1), mj1(ind2) 252 DO ji = 1, jpi 253 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea)