Changeset 12749
- Timestamp:
- 2020-04-15T15:49:42+02:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement
- Files:
-
- 9 edited
-
cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml (modified) (1 diff)
-
cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg (modified) (8 diffs)
-
cfgs/SHARED/field_def_nemo-oce.xml (modified) (1 diff)
-
cfgs/SHARED/namelist_ref (modified) (4 diffs)
-
cfgs/ref_cfgs.txt (modified) (1 diff)
-
src/ABL/abl.F90 (modified) (2 diffs)
-
src/ABL/ablmod.F90 (modified) (35 diffs)
-
src/ABL/par_abl.F90 (modified) (1 diff)
-
src/ABL/sbcabl.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml
r12063 r12749 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_r12563_ASINTER-06_ABL_improvement/cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg
r12489 r12749 22 22 cn_exp = "ORCA2" ! experience name 23 23 nn_it000 = 1 ! first time step 24 nn_itend = 16 ! last time step (std 5475)24 nn_itend = 5840 ! last time step (std 5840) 25 25 nn_date0 = 20130101 ! date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 26 26 nn_write = 4 ! frequency of write in the output file (modulo referenced to nn_it000) … … 31 31 &namdom ! time and space domain 32 32 !----------------------------------------------------------------------- 33 rn_Dt = 5400. ! time step for the dynamics and tracer33 rn_Dt = 5400. ! time step for the dynamics and tracer 34 34 / 35 35 !----------------------------------------------------------------------- … … 56 56 sn_sal = 'data_1m_salinity_nomask' , -1. ,'vosaline', .true. , .true. , 'yearly' , '' , '' , '' 57 57 / 58 58 59 !!====================================================================== 59 60 !! *** Surface Boundary Condition namelists *** !! … … 69 70 !! namsbc_rnf river runoffs (ln_rnf =T) 70 71 !! namsbc_apr Atmospheric Pressure (ln_apr_dyn =T) 71 !! namsbc_isf ice shelf melting/freezing (ln_isfcav =T : read (ln_read_cfg=T) or set or usr_def_zgr )72 !! namsbc_iscpl coupling option between land ice model and ocean (ln_isfcav =T)73 72 !! namsbc_wave external fields from wave model (ln_wave =T) 74 73 !! namberg iceberg floats (ln_icebergs=T) … … 79 78 !----------------------------------------------------------------------- 80 79 nn_fsbc = 1 ! frequency of SBC module call 81 ! (also = the frequency ofsea-ice & iceberg model call)82 ! Type of air-sea fluxes 80 ! ! (control sea-ice & iceberg model call) 81 ! Type of air-sea fluxes 83 82 ln_blk = .false. ! Bulk formulation (T => fill namsbc_blk ) 84 83 ln_abl = .true. ! ABL formulation (T => fill namsbc_abl ) 85 84 ! Sea-ice : 86 nn_ice = 2 ! =2 or 3 automatically for SI3 or CICE ("key_si3" or "key_cice") 87 ! except in AGRIF zoom where it has to be specified 88 ! Misc. options of sbc : 85 nn_ice = 2 ! =0 no ice boundary condition 86 ! ! =1 use observed ice-cover ( => fill namsbc_iif ) 87 ! ! =2 or 3 for SI3 and CICE, respectively 88 ln_ice_embd = .false. ! =T embedded sea-ice (pressure + mass and salt exchanges) 89 ! ! =F levitating ice (no pressure, mass and salt exchanges) 90 ! Misc. options of sbc : 89 91 ln_traqsr = .true. ! Light penetration in the ocean (T => fill namtra_qsr) 92 ln_dm2dc = .true. ! daily mean to diurnal cycle on short wave 90 93 ln_ssr = .true. ! Sea Surface Restoring on T and/or S (T => fill namsbc_ssr) 91 ln_dm2dc = .true. ! daily mean to diurnal cycle on short wave 94 nn_fwb = 2 ! FreshWater Budget: =0 unchecked 95 ! ! =1 global mean of e-p-r set to zero at each time step 96 ! ! =2 annual global mean of e-p-r set to zero 92 97 ln_rnf = .true. ! runoffs (T => fill namsbc_rnf) 93 nn_fwb = 2 ! FreshWater Budget:94 ! ! =2 annual global mean of e-p-r set to zero95 98 ln_wave = .false. ! Activate coupling with wave (T => fill namsbc_wave) 96 99 ln_cdgw = .false. ! Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave) 97 ln_sdw = .false. ! Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 100 ln_sdw = .false. ! Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 98 101 nn_sdrift = 0 ! Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift 99 102 ! ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] … … 111 114 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 112 115 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013) 113 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31) 114 rn_zqt = 10. ! Air temperature & humidity reference height (m) 115 rn_zu = 10. ! Wind vector reference height (m) 116 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 45r1) 116 117 ! 117 ! Skin is ONLY available in ECMWF and COARE algorithms: 118 ln_skin_cs = .false. ! use the cool-skin parameterization => set nn_fsbc=1 and ln_dm2dc=.true.! 119 ln_skin_wl = .false. ! use the warm-layer " => set nn_fsbc=1 and ln_dm2dc=.true.! 120 ! 121 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true. 122 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true. 123 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true. 118 ln_humi_sph = .true. ! humidity "sn_humi" is specific humidity [kg/kg] 119 ln_tpot = .false. !!GS: compute potential temperature or not 124 120 ! 125 cn_dir = './'! root directory for the bulk data location121 cn_dir = './' ! root directory for the bulk data location 126 122 !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! 127 123 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! … … 146 142 cn_dir = './' ! root directory for the location of the ABL grid file 147 143 cn_dom = 'dom_cfg_abl_L25Z10.nc' 144 145 ln_rstart_abl = .false. 148 146 ln_hpgls_frc = .true. 149 147 ln_geos_winds = .false. 150 148 nn_dyn_restore = 1 151 rn_ldyn_min = 7.5 ! magnitude of the nudging on ABL dynamics at the bottom of the ABL [hour] 149 150 rn_ldyn_min = 4.5 ! magnitude of the nudging on ABL dynamics at the bottom of the ABL [hour] 152 151 rn_ldyn_max = 1.5 ! magnitude of the nudging on ABL dynamics at the top of the ABL [hour] 153 rn_ltra_min = 7.5 ! magnitude of the nudging on ABL tracers at the bottom of the ABL [hour]152 rn_ltra_min = 4.5 ! magnitude of the nudging on ABL tracers at the bottom of the ABL [hour] 154 153 rn_ltra_max = 1.5 ! magnitude of the nudging on ABL tracers at the top of the ABL [hour] 155 154 / … … 201 200 &namberg ! iceberg parameters (default: OFF) 202 201 !----------------------------------------------------------------------- 202 ln_icebergs = .true. ! activate iceberg floats (force =F with "key_agrif") 203 204 cn_dir = './' ! root directory for the location of drag coefficient files 205 !______!___________!___________________!______________!______________!_________!___________!__________!__________!_______________! 206 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 207 ! ! ! (if <0 months) ! name ! (logical) ! (T/F ) ! 'monthly' ! filename ! pairing ! filename ! 208 sn_icb = 'calving', -1. , 'calving' , .true. , .true. , 'yearly' , '' , '' , '' 203 209 / 204 210 !!====================================================================== -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/SHARED/field_def_nemo-oce.xml
r12377 r12749 440 440 <field id="t_dta" long_name="DTA potential temperature" standard_name="dta_theta" unit="K" /> 441 441 <field id="q_dta" long_name="DTA specific humidity" standard_name="dta_qspe" unit="kg/kg" /> 442 <field id="coeft" long_name="ABL nudging coefficient" standard_name="coeft" unit="" /> 442 <field id="u_geo" long_name="GEO i-horizontal velocity" standard_name="geo_x_velocity" unit="m/s" /> 443 <field id="v_geo" long_name="GEO j-horizontal velocity" standard_name="geo_y_velocity" unit="m/s" /> 443 444 <field id="tke_abl" long_name="ABL turbulent kinetic energy" standard_name="abl_tke" unit="m2/s2" /> 444 445 <field id="avm_abl" long_name="ABL turbulent viscosity" standard_name="abl_avm" unit="m2/s" /> 445 446 <field id="avt_abl" long_name="ABL turbulent diffusivity" standard_name="abl_avt" unit="m2/s" /> 446 <field id="mxl_abl" long_name="ABL mixing length" standard_name="abl_mxl" unit="m" /> 447 <field id="mxlm_abl" long_name="ABL master mixing length" standard_name="abl_mxlm" unit="m" /> 448 <field id="mxld_abl" long_name="ABL dissipative mixing length" standard_name="abl_mxld" unit="m" /> 447 449 </field_group> 448 450 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/SHARED/namelist_ref
r12589 r12749 259 259 !----------------------------------------------------------------------- 260 260 ! ! bulk algorithm : 261 ln_NCAR = . true.! "NCAR" algorithm (Large and Yeager 2008)261 ln_NCAR = .false. ! "NCAR" algorithm (Large and Yeager 2008) 262 262 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 263 263 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013) … … 271 271 rn_pfac = 1. ! multipl. factor for precipitation (total & snow) 272 272 rn_efac = 1. ! multipl. factor for evaporation (0. or 1.) 273 rn_vfac = 0. ! multipl. factor for ocean & ice velocity273 rn_vfac = 1. ! multipl. factor for ocean & ice velocity 274 274 ! ! used to calculate the wind stress 275 275 ! ! (0. => absolute or 1. => relative winds) … … 277 277 ln_skin_wl = .false. ! use the warm-layer parameterization 278 278 ! ! ==> only available in ECMWF and COARE algorithms 279 ln_humi_sph = . true.! humidity "sn_humi" is specific humidity [kg/kg]279 ln_humi_sph = .false. ! humidity "sn_humi" is specific humidity [kg/kg] 280 280 ln_humi_dpt = .false. ! humidity "sn_humi" is dew-point temperature [K] 281 281 ln_humi_rlh = .false. ! humidity "sn_humi" is relative humidity [%] 282 ln_tpot = .true. !!GS: compute potential temperature or not 282 283 ! 283 284 cn_dir = './' ! root directory for the bulk data location … … 315 316 ! = 1 equatorial restoring 316 317 ! = 2 global restoring 317 rn_ldyn_min = 4.5 ! magnitude of the nudging on ABL dynamics at the bottom of the ABL [hour]318 rn_ldyn_max = 1.5 ! magnitude of the nudging on ABL dynamics at the top of the ABL [hour]319 rn_ltra_min = 4.5 ! magnitude of the nudging on ABL tracers at the bottom of the ABL [hour]320 rn_ltra_max = 1.5 ! magnitude of the nudging on ABL tracers at the top of the ABL [hour]318 rn_ldyn_min = 4.5 ! dynamics nudging magnitude inside the ABL [hour] (~3 rn_Dt) 319 rn_ldyn_max = 1.5 ! dynamics nudging magnitude above the ABL [hour] (~1 rn_Dt) 320 rn_ltra_min = 4.5 ! tracers nudging magnitude inside the ABL [hour] (~3 rn_Dt) 321 rn_ltra_max = 1.5 ! tracers nudging magnitude above the ABL [hour] (~1 rn_Dt) 321 322 nn_amxl = 0 ! mixing length: = 0 Deardorff 80 length-scale 322 323 ! = 1 length-scale based on the distance to the PBL height -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/ref_cfgs.txt
r12377 r12749 7 7 ORCA2_OFF_TRC OCE TOP OFF 8 8 ORCA2_SAS_ICE OCE ICE NST SAS 9 ORCA2_ICE_PISCES OCE TOP ICE NST 9 ORCA2_ICE_PISCES OCE TOP ICE NST ABL 10 10 ORCA2_ICE_ABL OCE ICE ABL 11 ORCA2_SAS_ICE_ABL OCE SAS ICE ABL12 ORCA2_ICE OCE ICE13 11 SPITZ12 OCE ICE 14 12 WED025 OCE ICE 15 eORCA025_ICE OCE ICE16 eORCA025_ICE_ABL OCE ICE ABL17 eORCA025_SAS_ICE_ABL OCE SAS ICE ABL -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/abl.F90
r12588 r12749 39 39 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: msk_abl 40 40 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: rest_eq 41 42 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: cft_abl43 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,: ) :: taux_abl44 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,: ) :: tauy_abl45 41 ! 46 42 INTEGER , PUBLIC :: nt_n, nt_a !: now / after indices (equal 1 or 2) … … 68 64 & mxld_abl(1:jpi,1:jpj,1:jpka ), & 69 65 & mxlm_abl(1:jpi,1:jpj,1:jpka ), & 70 & cft_abl (1:jpi,1:jpj,1:jpka ), &71 66 & fft_abl (1:jpi,1:jpj ), & 72 67 & pblh (1:jpi,1:jpj ), & 73 & taux_abl(1:jpi,1:jpj ), &74 & tauy_abl(1:jpi,1:jpj ), &75 68 & msk_abl (1:jpi,1:jpj ), & 76 69 & rest_eq (1:jpi,1:jpj ), & -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90
r12592 r12749 38 38 39 39 !=================================================================================================== 40 SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &! in40 SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, & ! in 41 41 & pu_dta, pv_dta, pt_dta, pq_dta, & 42 42 & pslp_dta, pgu_dta, pgv_dta, & 43 & pcd_du, psen, pevp, & 43 & pcd_du, psen, pevp, & ! in/out 44 44 & pwndm, ptaui, ptauj, ptaum & 45 45 #if defined key_si3 … … 112 112 REAL(wp) :: zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2 113 113 REAL(wp) :: zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce 114 LOGICAL :: SemiImp_Cor = . FALSE.114 LOGICAL :: SemiImp_Cor = .TRUE. 115 115 ! 116 116 !!--------------------------------------------------------------------- … … 135 135 ustar2(ji,jj) = zzoce 136 136 #endif 137 zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT( Cdn_oce(ji,jj)) ) !<-- recover the value of z0 from Cdn_oce137 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 138 138 END_2D 139 139 ! … … 162 162 DO ji = 1, jpi ! vector opt. 163 163 ! Neumann at the bottom 164 z_elem_a( ji, 2 ) = 0._wp165 z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2)164 z_elem_a( ji, 2 ) = 0._wp 165 z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 ) 166 166 ! Homogeneous Neumann at the top 167 167 z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 168 168 z_elem_c( ji, jpka ) = 0._wp 169 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, 169 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 170 170 END DO 171 171 … … 173 173 174 174 DO jk = 3, jpkam1 175 DO ji = 2, jpim1176 !DO ji = 1,jpi !!GS: to be checked177 tq_abl ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side175 !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 178 178 END DO 179 179 END DO … … 187 187 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj) 188 188 #endif 189 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2) + rDt_abl * zztmp1190 tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2189 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 191 191 tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra ) 192 192 END DO … … 199 199 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj) 200 200 #endif 201 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2) + rDt_abl * zztmp1202 tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2201 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 203 203 tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra ) 204 204 END DO … … 208 208 !! ---------------------------------------------------------- 209 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)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 213 END DO 214 214 215 215 DO jk = 3, jpka 216 216 DO ji = 1,jpi 217 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )217 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF( ji, jk-1 ) ) 218 218 zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) 219 219 tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk ,nt_a,jtra) & 220 & - z_elem_a(ji, jk) *tq_abl(ji,jj,jk-1,nt_a,jtra) )220 & - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) ) 221 221 END DO 222 222 END DO … … 240 240 IF( SemiImp_Cor ) THEN 241 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) 254 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 ELSE 279 280 !------------- 281 DO jk = 2, jpka ! outer loop 282 !------------- 283 ! 284 IF( MOD( kt, 2 ) == 0 ) then 285 ! Advance u_abl & v_abl to time n+1 286 DO jj = 1, jpj 287 DO ji = 1, jpi 288 zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_n ) - pgv_dta(ji ,jj ,jk) ) 289 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff 290 zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_a ) - pgu_dta(ji ,jj ,jk) ) 291 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff ) 292 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * u_abl( ji, jj, jk, nt_a ) 293 END DO 294 END DO 295 ! 296 ELSE 297 ! 298 DO jj = 1, jpj 299 DO ji = 1, jpi 300 zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_n ) - pgu_dta(ji ,jj ,jk) ) 301 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff 302 zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_a ) - pgv_dta(ji ,jj ,jk) ) 303 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff ) 304 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * v_abl( ji, jj, jk, nt_a ) 305 END DO 306 END DO 307 ! 308 ENDIF 309 ! 310 !------------- 311 END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 312 !------------- 313 314 ENDIF 315 316 !------------- 317 ! 318 IF( ln_hpgls_frc ) THEN 319 DO jj = 1, jpj ! outer loop 320 DO jk = 1, jpka 321 DO ji = 1, jpi 322 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 323 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 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) 254 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 324 285 ENDDO 325 286 ENDDO 326 ENDDO 327 END IF 328 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 329 327 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 330 328 ! ! 4 *** Advance u,v to time n+1 … … 338 336 DO jk = 3, jpkam1 339 337 DO ji = 1, jpi 340 z_elem_a( ji, 341 z_elem_c( ji, 342 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )! diagonal338 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 343 341 END DO 344 342 END DO … … 346 344 DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi) 347 345 !++ Surface boundary condition 348 z_elem_a( ji, 2) = 0._wp349 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2)346 z_elem_a( ji, 2 ) = 0._wp 347 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 350 348 ! 351 349 zztmp1 = pcd_du(ji, jj) … … 359 357 u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2 360 358 361 IF( ln_topbc_neumann ) THEN 362 !++ Top Neumann B.C. 363 z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 364 z_elem_c( ji, jpka ) = 0._wp 365 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 366 !u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl ( ji, jj, jpka, nt_a ) 367 ELSE 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 368 367 !++ Top Dirichlet B.C. 369 368 z_elem_a( ji, jpka ) = 0._wp … … 371 370 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 372 371 u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 373 ENDIF372 !ENDIF 374 373 375 374 END DO … … 377 376 !! Matrix inversion 378 377 !! ---------------------------------------------------------- 379 DO ji = 2, jpi 378 !DO ji = 2, jpi 379 DO ji = 1, jpi !!GS: TBI 380 380 zcff = 1._wp / z_elem_b( ji, 2 ) 381 381 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) … … 410 410 DO jk = 3, jpkam1 411 411 DO ji = 1, jpi 412 z_elem_a( ji, 413 z_elem_c( ji, 414 z_elem_b( ji, 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 415 END DO 416 416 END DO … … 418 418 DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 419 419 !++ Surface boundary condition 420 z_elem_a( ji, 2) = 0._wp421 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2)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 422 ! 423 423 zztmp1 = pcd_du(ji, jj) … … 431 431 v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2 432 432 433 IF( ln_topbc_neumann ) THEN 434 !++ Top Neumann B.C. 435 z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 436 z_elem_c( ji, jpka ) = 0._wp 437 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 438 !v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl ( ji, jj, jpka, nt_a ) 439 ELSE 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 440 441 !++ Top Dirichlet B.C. 441 442 z_elem_a( ji, jpka ) = 0._wp … … 443 444 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 444 445 v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 445 ENDIF446 !ENDIF 446 447 447 448 END DO … … 528 529 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 529 530 ! 530 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1., v_abl(:,:,:,nt_a) , 'T', -1. , kfillmode = jpfillnothing)531 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1., v_abl(:,:,:,nt_a) , 'T', -1. ) 531 532 CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T', 1., kfillmode = jpfillnothing ) ! ++++ this should not be needed... 532 533 ! 533 534 #if defined key_iomput 534 ! first ABL level 535 ! 2D & first ABL level 536 IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh", pblh(:,: ) ) 535 537 IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl", u_abl(:,:,2,nt_a ) ) 536 538 IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl", v_abl(:,:,2,nt_a ) ) … … 541 543 IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta", pt_dta(:,:,2 ) ) 542 544 IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta", pq_dta(:,:,2 ) ) 543 ! all ABL levels 544 IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) ) 545 IF ( iom_use("v_abl" ) ) CALL iom_put ( "v_abl" , v_abl(:,:,2:jpka,nt_a ) ) 546 IF ( iom_use("t_abl" ) ) CALL iom_put ( "t_abl" , tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 547 IF ( iom_use("q_abl" ) ) CALL iom_put ( "q_abl" , tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 548 IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a ) ) 549 IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka ) ) 550 IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka ) ) 551 IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxlm_abl(:,:,2:jpka ) ) 552 IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh" , pblh(:,: ) ) 553 ! 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 554 565 IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta", pu_dta(:,:,2:jpka) ) 555 566 IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta", pv_dta(:,:,2:jpka) ) … … 557 568 IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta", pq_dta(:,:,2:jpka) ) 558 569 IF( ln_geos_winds ) THEN 559 IF ( iom_use("u z1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) )560 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) ) 561 572 END IF 562 573 IF( ln_hpgls_frc ) THEN 563 IF ( iom_use("u z1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp))564 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) ) 565 576 END IF 566 577 #endif … … 570 581 ! 571 582 DO_2D_11_11 572 ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta )573 zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa )583 ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta ) 584 zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa ) 574 585 zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 575 586 psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) !GS: negative sign to respect aerobulk convention … … 587 598 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 588 599 DO_2D_11_11 589 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 590 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 591 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 592 pwndm (ji,jj) = zcff 593 ptaum (ji,jj) = zztmp * zcff 594 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 595 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 596 taux_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_i(ji,jj) 597 tauy_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_j(ji,jj) 600 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 601 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 602 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 603 604 pwndm (ji,jj) = zcff 605 ptaum (ji,jj) = zztmp * zcff 606 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 607 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 598 608 END_2D 599 609 ! ... utau, vtau at U- and V_points, resp. … … 620 630 621 631 #if defined key_si3 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 632 ! ------------------------------------------------------------ ! 633 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 634 ! ------------------------------------------------------------ ! 635 DO_2D_00_00 636 637 zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 638 zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 639 640 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 641 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 642 & * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 643 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 644 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 645 & * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 646 END_2D 647 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 648 ! 649 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' & 650 & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' ) 641 651 #endif 642 652 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 688 698 REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv 689 699 REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2 690 REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 700 REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 ! zbuoy for BL89 691 701 REAL(wp) :: zwndi, zwndj 692 702 REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2 … … 733 743 ! Terms for the tridiagonal problem 734 744 DO jk = 2, jpkam1 735 DO ji = 1, jpi736 zshear = zsh2( ji, jk )! zsh2 is already multiplied by Avm_abl at this point737 zsh2(ji,jk) 738 zbuoy 739 740 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-diagonal741 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-diagonal745 DO ji = 1, jpi 746 zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point 747 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 748 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 749 750 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 751 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 742 752 IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE 743 z_elem_b( ji, jk )= e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &744 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)! diagonal745 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) ) 753 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 754 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) ! diagonal 755 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) ) ! right-hand-side 746 756 ELSE 747 z_elem_b( ji, jk )= e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &748 & 749 & 750 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear )! right-hand-side757 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 758 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) & ! diagonal 759 & - e3w_abl(jk) * rDt_abl * zbuoy 760 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side 751 761 END IF 752 762 END DO … … 754 764 755 765 DO ji = 1,jpi ! vector opt. 756 zesrf = MAX( rn_esfc * ustar2(ji,jj), tke_min ) 757 zetop = tke_min 758 759 z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp 760 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 ) 761 z_elem_c ( ji, jpka ) = 0._wp 762 z_elem_b ( ji, jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 763 764 tke_abl ( ji, jj, 1 , nt_a ) = zesrf 765 tke_abl ( ji, jj, jpka , nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 766 767 zbn2(ji,jj, 1) = zbn2( ji,jj, 2) 768 zsh2(ji, 1) = zsh2( ji, 2) 769 zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 770 zsh2(ji, jpka) = zsh2( ji , jpkam1) 766 zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min ) 767 zetop = tke_min 768 769 z_elem_a ( ji, 1 ) = 0._wp 770 z_elem_c ( ji, 1 ) = 0._wp 771 z_elem_b ( ji, 1 ) = 1._wp 772 tke_abl ( ji, jj, 1, nt_a ) = zesrf 773 774 !++ Top Neumann B.C. 775 !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 ) 776 !z_elem_c ( ji, jpka ) = 0._wp 777 !z_elem_b ( ji, jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 778 !tke_abl ( ji, jj, jpka, nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 779 780 !++ Top Dirichlet B.C. 781 z_elem_a ( ji, jpka ) = 0._wp 782 z_elem_c ( ji, jpka ) = 0._wp 783 z_elem_b ( ji, jpka ) = 1._wp 784 tke_abl ( ji, jj, jpka, nt_a ) = zetop 785 786 zbn2 ( ji, jj, 1 ) = zbn2 ( ji, jj, 2 ) 787 zsh2 ( ji, 1 ) = zsh2 ( ji, 2 ) 788 zbn2 ( ji, jj, jpka ) = zbn2 ( ji, jj, jpkam1 ) 789 zsh2 ( ji, jpka ) = zsh2 ( ji , jpkam1 ) 771 790 END DO 772 791 !! … … 844 863 ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 845 864 IF(ln_smth_pblh) THEN 846 CALL lbc_lnk( 'ablmod', pblh, 'T', 1. , kfillmode = jpfillnothing)865 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing) 847 866 CALL smooth_pblh( pblh, msk_abl ) 848 CALL lbc_lnk( 'ablmod', pblh, 'T', 1. , kfillmode = jpfillnothing)867 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing) 849 868 ENDIF 850 869 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 860 879 ! 861 880 DO ji = 1, jpi 862 mxld_abl 863 mxld_abl 864 mxlm_abl 865 mxlm_abl 866 zldw 867 zlup 881 mxld_abl( ji, jj, 1 ) = mxl_min 882 mxld_abl( ji, jj, jpka ) = mxl_min 883 mxlm_abl( ji, jj, 1 ) = mxl_min 884 mxlm_abl( ji, jj, jpka ) = mxl_min 885 zldw ( ji, 1 ) = zrough(ji,jj) * rn_Lsfc 886 zlup ( ji, jpka ) = mxl_min 868 887 END DO 869 888 ! … … 898 917 DO jk = 1, jpka 899 918 DO ji = 1, jpi 919 ! zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 900 920 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 901 ! zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp)902 921 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 903 !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) )904 922 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 905 923 END DO … … 971 989 # undef zlup 972 990 # undef zldw 973 !!991 ! 974 992 CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale 975 993 ! … … 1007 1025 DO jkup=jk+1,jpka-1 1008 1026 DO ji = 1, jpi 1009 zbuoy1 1010 zbuoy2 1027 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1028 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 1011 1029 zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 1012 1030 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & … … 1038 1056 DO jkdwn=jk-1,1,-1 1039 1057 DO ji = 1, jpi 1040 zbuoy1 1041 zbuoy2 1058 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1059 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 1042 1060 zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) & 1043 1061 & * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & … … 1166 1184 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1167 1185 zldw(ji,jk) = zcff 1168 1186 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1169 1187 ln_foundl(ji) = .true. 1170 1188 END IF … … 1186 1204 # undef zlup 1187 1205 # undef zldw 1188 !1189 1206 ! 1190 1207 ! -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/par_abl.F90
r12588 r12749 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. !: smoothing of atmospheric PBL height30 !LOGICAL , PUBLIC :: ln_topbc_neumann = .FALSE. !: idealised testcases only 31 31 32 32 LOGICAL , PUBLIC :: ln_rstart_abl !: (de)activate abl restart -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90
r12592 r12749 168 168 rn_Esfc = 1._wp / SQRT(rn_cm*rn_ceps) 169 169 rn_Lsfc = vkarmn * SQRT(SQRT(rn_cm*rn_ceps)) / rn_cm 170 170 171 171 IF(lwp) THEN 172 172 WRITE(numout,*) … … 175 175 IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 176 176 IF(nn_amxl==1) WRITE(numout,*) 'Modified Deardorff 80 length-scale ' 177 IF(nn_amxl== 1) WRITE(numout,*) 'Bougeault and Lacarrere length-scale '177 IF(nn_amxl==2) WRITE(numout,*) 'Bougeault and Lacarrere length-scale ' 178 178 IF(nn_amxl==3) WRITE(numout,*) 'Rodier et al. length-scale ' 179 179 WRITE(numout,*) ' Minimum value of atmospheric TKE = ',tke_min,' m^2 s^-2'
Note: See TracChangeset
for help on using the changeset viewer.
