Changeset 12154 for NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019
- Timestamp:
- 2019-12-10T15:44:23+01:00 (16 months ago)
- Location:
- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019
- Files:
-
- 2 deleted
- 44 edited
- 10 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg
r11536 r12154 93 93 !----------------------------------------------------------------------- 94 94 ! ! bulk algorithm : 95 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 96 95 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 96 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 97 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013) 98 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31) 99 ! 100 rn_zqt = 10. ! Air temperature & humidity reference height (m) 101 rn_zu = 10. ! Wind vector reference height (m) 102 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012) 103 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015) 104 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow) 105 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.) 106 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to 107 ! ! calculate the wind stress (0.=absolute or 1.=relative winds) 108 ln_skin_cs = .false. ! use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB 109 ln_skin_wl = .false. ! use the warm-layer " " " 110 ! 111 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true. 112 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true. 113 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true. 114 ! 97 115 cn_dir = './' ! root directory for the bulk data location 98 116 !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! … … 108 126 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 109 127 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 110 sn_tdif = 'taudif_core' , 24. , 'taudif' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , ''111 128 / 112 129 !----------------------------------------------------------------------- -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
r11536 r12154 89 89 !----------------------------------------------------------------------- 90 90 ! ! bulk algorithm : 91 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 92 91 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 92 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 93 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013) 94 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31) 95 ! 96 rn_zqt = 10. ! Air temperature & humidity reference height (m) 97 rn_zu = 10. ! Wind vector reference height (m) 98 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012) 99 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015) 100 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow) 101 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.) 102 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to 103 ! ! calculate the wind stress (0.=absolute or 1.=relative winds) 104 ln_skin_cs = .false. ! use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB 105 ln_skin_wl = .false. ! use the warm-layer " " " 106 ! 107 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true. 108 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true. 109 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true. 110 ! 93 111 cn_dir = './' ! root directory for the bulk data location 94 112 !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! … … 104 122 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core2_nordic1_bilin.nc' , '' , '' 105 123 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core2_nordic1_bilin.nc' , '' , '' 106 sn_tdif = 'taudif_core' , 24. , 'taudif' , .false. , .true. , 'yearly' , 'weights_core2_nordic1_bilin.nc' , '' , ''107 124 108 125 / -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/AGRIF_DEMO/EXPREF/2_namelist_ice_ref
r9575 r12154 1 link 1_namelist_ice_ref1 link ../../SHARED/namelist_ice_ref -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/AGRIF_DEMO/EXPREF/2_namelist_ref
r9464 r12154 1 link 1_namelist_ref1 link ../../SHARED/namelist_ref -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
r11536 r12154 104 104 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core2_nordic2_bilin.nc' , '' , '' 105 105 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core2_nordic2_bilin.nc' , '' , '' 106 sn_tdif = 'taudif_core' , 24. , 'taudif' , .false. , .true. , 'yearly' , 'weights_core2_nordic2_bilin.nc' , '' , ''107 106 108 107 / -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
r11536 r12154 93 93 !----------------------------------------------------------------------- 94 94 ! ! bulk algorithm : 95 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 96 95 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 96 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 97 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013) 98 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31) 99 ! 100 rn_zqt = 10. ! Air temperature & humidity reference height (m) 101 rn_zu = 10. ! Wind vector reference height (m) 102 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012) 103 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015) 104 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow) 105 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.) 106 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to 107 ! ! calculate the wind stress (0.=absolute or 1.=relative winds) 108 ln_skin_cs = .false. ! use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB 109 ln_skin_wl = .false. ! use the warm-layer " " " 110 ! 111 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true. 112 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true. 113 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true. 114 ! 97 115 cn_dir = './' ! root directory for the bulk data location 98 116 !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! … … 108 126 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 109 127 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 110 sn_tdif = 'taudif_core' , 24. , 'taudif' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , ''111 128 / 112 129 !----------------------------------------------------------------------- -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/C1D_PAPA/EXPREF/namelist_cfg
r11536 r12154 159 159 sn_snow = 'forcing_C1D_PAPA' , 3. , 'sososnow', .false. , .false. , 'yearly' , '' , '' , '' 160 160 sn_slp = 'forcing_C1D_PAPA' , 3. , 'somslpre', .true. , .false. , 'yearly' , '' , '' , '' 161 sn_tdif = 'forcing_C1D_PAPA' , 24. , 'taudif' , .false. , .false. , 'yearly' , '' , '' , ''162 161 163 162 / -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml
r11306 r12154 9 9 --> 10 10 11 <file_definition type="one_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq=" 6h" min_digits="4">11 <file_definition type="one_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="5d" min_digits="4"> 12 12 13 <file_group id=" 6h" output_freq="6h" output_level="10" enabled=".TRUE."> <!-- 6hfiles -->13 <file_group id="5d" output_freq="5d" output_level="10" enabled=".TRUE."> <!-- 5d files --> 14 14 <file id="file11" name_suffix="_grid_T" description="ocean T grid variables" > 15 15 <field field_ref="e3t" /> -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg
r11945 r12154 89 89 ln_traqsr = .true. ! Light penetration in the ocean (T => fill namtra_qsr) 90 90 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 91 92 ln_rnf = .true. ! runoffs (T => fill namsbc_rnf) 92 93 nn_fwb = 2 ! FreshWater Budget: … … 107 108 !----------------------------------------------------------------------- 108 109 ! ! bulk algorithm : 109 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 110 110 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 111 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 112 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 ! 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. 124 ! 111 125 cn_dir = './' ! root directory for the bulk data location 112 126 !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! … … 404 418 !! namdiu Cool skin and warm layer models (default: OFF) 405 419 !! namdiu Cool skin and warm layer models (default: OFF) 406 !! namflo float parameters ( "key_float")407 !! nam_diaharm Harmonic analysis of tidal constituents ( "key_diaharm")408 !! nam dct transports through some sections ("key_diadct")420 !! namflo float parameters (default: OFF) 421 !! nam_diaharm Harmonic analysis of tidal constituents (default: OFF) 422 !! nam_diadct transports through some sections (default: OFF) 409 423 !! nam_diatmb Top Middle Bottom Output (default: OFF) 410 424 !! nam_dia25h 25h Mean Output (default: OFF) -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
r11536 r12154 118 118 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 119 119 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 120 sn_tdif = 'taudif_core' , 24. , 'taudif' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , ''121 120 / 122 121 !----------------------------------------------------------------------- -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg
r11536 r12154 64 64 &namsbc_blk ! namsbc_blk generic Bulk formula (ln_blk =T) 65 65 !----------------------------------------------------------------------- 66 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 67 66 ! ! bulk algorithm : 67 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 68 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 69 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013) 70 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31) 71 ! 72 rn_zqt = 10. ! Air temperature & humidity reference height (m) 73 rn_zu = 10. ! Wind vector reference height (m) 74 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012) 75 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015) 76 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow) 77 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.) 78 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to 79 ! ! calculate the wind stress (0.=absolute or 1.=relative winds) 80 ln_skin_cs = .false. ! use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB 81 ln_skin_wl = .false. ! use the warm-layer " " " 82 ! 83 ln_humi_sph = .true. ! humidity specified below in "sn_humi" is specific humidity [kg/kg] if .true. 84 ln_humi_dpt = .false. ! humidity specified below in "sn_humi" is dew-point temperature [K] if .true. 85 ln_humi_rlh = .false. ! humidity specified below in "sn_humi" is relative humidity [%] if .true. 86 ! 68 87 cn_dir = './' ! root directory for the bulk data location 69 88 !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! … … 79 98 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 80 99 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 81 sn_tdif = 'taudif_core' , 24. , 'taudif' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , ''82 100 / 83 101 !----------------------------------------------------------------------- -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/SHARED/field_def_nemo-oce.xml
r12109 r12154 41 41 <field id="ahmt_3d" long_name=" 3D t-eddy viscosity coefficient" unit="m2/s or m4/s" grid_ref="grid_T_3D"/> 42 42 43 <field id="sst" long_name="sea surface temperature" standard_name="sea_surface_temperature" unit="degC" /> 43 <field id="sst" long_name="Bulk sea surface temperature" standard_name="bulk_sea_surface_temperature" unit="degC" /> 44 <field id="t_skin" long_name="Skin temperature aka SSST" standard_name="skin_temperature" unit="degC" /> 44 45 <field id="sst2" long_name="square of sea surface temperature" standard_name="square_of_sea_surface_temperature" unit="degC2" > sst * sst </field > 45 46 <field id="sstmax" long_name="max of sea surface temperature" field_ref="sst" operation="maximum" /> … … 300 301 301 302 <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> 303 <field id="rho_air" long_name="Air density at 10m above sea surface" standard_name="rho_air_10m" unit="kg/m3" /> 304 <field id="dt_skin" long_name="SSST-SST temperature difference" standard_name="SSST-SST" unit="K" /> 302 305 <field id="qlw_oce" long_name="Longwave Downward Heat Flux over open ocean" standard_name="surface_net_downward_longwave_flux" unit="W/m2" /> 303 306 <field id="qsb_oce" long_name="Sensible Downward Heat Flux over open ocean" standard_name="surface_downward_sensible_heat_flux" unit="W/m2" /> 304 307 <field id="qla_oce" long_name="Latent Downward Heat Flux over open ocean" standard_name="surface_downward_latent_heat_flux" unit="W/m2" /> 308 <field id="evap_oce" long_name="Evaporation over open ocean" standard_name="evaporation" unit="kg/m2/s" /> 305 309 <field id="qt_oce" long_name="total flux at ocean surface" standard_name="surface_downward_heat_flux_in_sea_water" unit="W/m2" /> 306 310 <field id="qsr_oce" long_name="solar heat flux at ocean surface" standard_name="net_downward_shortwave_flux_at_sea_water_surface" unit="W/m2" /> … … 362 366 </field_group> 363 367 364 <!-- scalar variables --> 365 <field_group id="SBC_0D" grid_ref="grid_1point" > 368 369 </field_group> <!-- SBC --> 370 371 <!-- ABL --> 372 <field_group id="ABL" > <!-- time step automaticaly defined based on nn_fsbc --> 373 374 <!-- variables available with ABL on atmospheric T grid--> 375 <field_group id="grid_ABL3D" grid_ref="grid_TA_3D" > 376 <field id="u_abl" long_name="ABL i-horizontal velocity" standard_name="abl_x_velocity" unit="m/s" /> 377 <field id="v_abl" long_name="ABL j-horizontal velocity" standard_name="abl_y_velocity" unit="m/s" /> 378 <field id="t_abl" long_name="ABL potential temperature" standard_name="abl_theta" unit="K" /> 379 <field id="q_abl" long_name="ABL specific humidity" standard_name="abl_qspe" unit="kg/kg" /> 380 <!-- debug (to be removed) --> 381 <field id="u_dta" long_name="DTA i-horizontal velocity" standard_name="dta_x_velocity" unit="m/s" /> 382 <field id="v_dta" long_name="DTA j-horizontal velocity" standard_name="dta_y_velocity" unit="m/s" /> 383 <field id="t_dta" long_name="DTA potential temperature" standard_name="dta_theta" unit="K" /> 384 <field id="q_dta" long_name="DTA specific humidity" standard_name="dta_qspe" unit="kg/kg" /> 385 <field id="coeft" long_name="ABL nudging coefficient" standard_name="coeft" unit="" /> 386 <field id="tke_abl" long_name="ABL turbulent kinetic energy" standard_name="abl_tke" unit="m2/s2" /> 387 <field id="avm_abl" long_name="ABL turbulent viscosity" standard_name="abl_avm" unit="m2/s" /> 388 <field id="avt_abl" long_name="ABL turbulent diffusivity" standard_name="abl_avt" unit="m2/s" /> 389 <field id="mxl_abl" long_name="ABL mixing length" standard_name="abl_mxl" unit="m" /> 366 390 </field_group> 367 391 368 </field_group> <!-- SBC --> 369 392 <field_group id="grid_ABL2D" grid_ref="grid_TA_2D" > 393 <field id="pblh" long_name="ABL height" standard_name="abl_height" unit="m" /> 394 <field id="uz1_abl" long_name="ABL i-horizontal velocity" standard_name="abl_x_velocity" unit="m/s" /> 395 <field id="vz1_abl" long_name="ABL j-horizontal velocity" standard_name="abl_y_velocity" unit="m/s" /> 396 <field id="uvz1_abl" long_name="ABL wind speed module" standard_name="abl_wind_speed" unit="m/s" > sqrt( uz1_abl^2 + vz1_abl^2 ) </field> 397 <field id="tz1_abl" long_name="ABL potential temperature" standard_name="abl_theta" unit="K" /> 398 <field id="qz1_abl" long_name="ABL specific humidity" standard_name="abl_qspe" unit="kg/kg" /> 399 <field id="uz1_dta" long_name="DTA i-horizontal velocity" standard_name="dta_x_velocity" unit="m/s" /> 400 <field id="vz1_dta" long_name="DTA j-horizontal velocity" standard_name="dta_y_velocity" unit="m/s" /> 401 <field id="uvz1_dta" long_name="DTA wind speed module" standard_name="dta_wind_speed" unit="m/s" > sqrt( uz1_dta^2 + vz1_dta^2 ) </field> 402 <field id="tz1_dta" long_name="DTA potential temperature" standard_name="dta_theta" unit="K" /> 403 <field id="qz1_dta" long_name="DTA specific humidity" standard_name="dta_qspe" unit="kg/kg" /> 404 <!-- debug (to be removed) --> 405 <field id="uz1_geo" long_name="GEO i-horizontal velocity" standard_name="geo_x_velocity" unit="m/s" /> 406 <field id="vz1_geo" long_name="GEO j-horizontal velocity" standard_name="geo_y_velocity" unit="m/s" /> 407 <field id="uvz1_geo" long_name="GEO wind speed module" standard_name="geo_wind_speed" unit="m/s" > sqrt( uz1_geo^2 + vz1_geo^2 ) </field> 408 </field_group> 409 410 </field_group> <!-- ABL --> 411 412 370 413 <!-- U grid --> 371 414 … … 391 434 <field id="uocet" long_name="ocean transport along i-axis times temperature (CRS)" unit="degC*m/s" grid_ref="grid_U_3D" /> 392 435 <field id="uoces" long_name="ocean transport along i-axis times salinity (CRS)" unit="1e-3*m/s" grid_ref="grid_U_3D" /> 436 <field id="ssuww" long_name="ocean surface wind work along i-axis" standard_name="surface_x_wind_work" unit="N/m*s" > utau * ssu </field> 393 437 394 438 <!-- u-eddy diffusivity coefficients (available if ln_traldf_OFF=F) --> … … 448 492 <field id="vocet" long_name="ocean transport along j-axis times temperature (CRS)" unit="degC*m/s" grid_ref="grid_V_3D" /> 449 493 <field id="voces" long_name="ocean transport along j-axis times salinity (CRS)" unit="1e-3*m/s" grid_ref="grid_V_3D" /> 494 <field id="ssvww" long_name="ocean surface wind work along j-axis" standard_name="surface_y_wind_work" unit="N/m*s" > vtau * ssv </field> 450 495 451 496 <!-- v-eddy diffusivity coefficients (available if ln_traldf_OFF=F) --> … … 589 634 590 635 591 <!-- variables available with key_float-->636 <!-- variables available with ln_floats --> 592 637 593 638 <field_group id="floatvar" grid_ref="grid_T_nfloat" operation="instant" > -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/SHARED/namelist_ref
r12113 r12154 5 5 !! namelists 2 - Surface boundary (namsbc, namsbc_flx, namsbc_blk, namsbc_cpl, 6 6 !! namsbc_sas, namtra_qsr, namsbc_rnf, 7 !! namsbc_isf, namsbc_iscpl, namsbc_apr, 7 !! namsbc_isf, namsbc_iscpl, namsbc_apr, 8 8 !! namsbc_ssr, namsbc_wave, namberg) 9 9 !! 3 - lateral boundary (namlbc, namagrif, nambdy, nambdy_tide) … … 65 65 ln_clobber = .true. ! clobber (overwrite) an existing file 66 66 nn_chunksz = 0 ! chunksize (bytes) for NetCDF file (works only with iom_nf90 routines) 67 ln_xios_read = . FALSE. ! use XIOS to read restart file (only for a single file restart)67 ln_xios_read = .false. ! use XIOS to read restart file (only for a single file restart) 68 68 nn_wxios = 0 ! use XIOS to write restart file 0 - no, 1 - single file output, 2 - multiple file output 69 69 / … … 88 88 cn_domcfg = "domain_cfg" ! domain configuration filename 89 89 ! 90 ln_closea = .false. ! T => keep closed seas (defined by closea_mask field) in the 90 ln_closea = .false. ! T => keep closed seas (defined by closea_mask field) in the 91 91 ! ! domain and apply special treatment of freshwater fluxes. 92 ! ! F => suppress closed seas (defined by closea_mask field) 92 ! ! F => suppress closed seas (defined by closea_mask field) 93 93 ! ! from the bathymetry at runtime. 94 94 ! ! If closea_mask field doesn't exist in the domain_cfg file … … 106 106 ln_tsd_init = .false. ! ocean initialisation 107 107 ln_tsd_dmp = .false. ! T-S restoring (see namtra_dmp) 108 108 109 109 cn_dir = './' ! root directory for the T-S data location 110 110 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! … … 195 195 nn_fsbc = 2 ! frequency of SBC module call 196 196 ! ! (control sea-ice & iceberg model call) 197 ! Type of air-sea fluxes 197 ! Type of air-sea fluxes 198 198 ln_usr = .false. ! user defined formulation (T => check usrdef_sbc) 199 199 ln_flx = .false. ! flux formulation (T => fill namsbc_flx ) 200 200 ln_blk = .false. ! Bulk formulation (T => fill namsbc_blk ) 201 ln_abl = .false. ! ABL formulation (T => fill namsbc_abl ) 201 202 ! ! Type of coupling (Ocean/Ice/Atmosphere) : 202 203 ln_cpl = .false. ! atmosphere coupled formulation ( requires key_oasis3 ) … … 205 206 ! ! =0 no opa-sas OASIS coupling: default single executable config. 206 207 ! ! =1 opa-sas OASIS coupling: multi executable config., OPA component 207 ! ! =2 opa-sas OASIS coupling: multi executable config., SAS component 208 ! ! =2 opa-sas OASIS coupling: multi executable config., SAS component 208 209 ! Sea-ice : 209 nn_ice = 0 ! =0 no ice boundary condition 210 nn_ice = 0 ! =0 no ice boundary condition 210 211 ! ! =1 use observed ice-cover ( => fill namsbc_iif ) 211 212 ! ! =2 or 3 automatically for SI3 or CICE ("key_si3" or "key_cice") … … 213 214 ln_ice_embd = .false. ! =T embedded sea-ice (pressure + mass and salt exchanges) 214 215 ! ! =F levitating ice (no pressure, mass and salt exchanges) 215 ! Misc. options of sbc : 216 ! Misc. options of sbc : 216 217 ln_traqsr = .false. ! Light penetration in the ocean (T => fill namtra_qsr) 217 218 ln_dm2dc = .false. ! daily mean to diurnal cycle on short wave … … 225 226 ln_wave = .false. ! Activate coupling with wave (T => fill namsbc_wave) 226 227 ln_cdgw = .false. ! Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave) 227 ln_sdw = .false. ! Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 228 ln_sdw = .false. ! Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 228 229 nn_sdrift = 0 ! Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift 229 230 ! ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] … … 250 251 / 251 252 !----------------------------------------------------------------------- 252 &namsbc_blk ! namsbc_blk generic Bulk formula 253 &namsbc_blk ! namsbc_blk generic Bulk formula (ln_blk =T) 253 254 !----------------------------------------------------------------------- 254 255 ! ! bulk algorithm : 255 ln_NCAR = .false.! "NCAR" algorithm (Large and Yeager 2008)256 ln_NCAR = .true. ! "NCAR" algorithm (Large and Yeager 2008) 256 257 ln_COARE_3p0 = .false. ! "COARE 3.0" algorithm (Fairall et al. 2003) 257 ln_COARE_3p 5 = .false. ! "COARE 3.5" algorithm (Edson et al. 2013)258 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 31)258 ln_COARE_3p6 = .false. ! "COARE 3.6" algorithm (Edson et al. 2013) 259 ln_ECMWF = .false. ! "ECMWF" algorithm (IFS cycle 45r1) 259 260 ! 260 rn_zqt = 10. ! Air temperature & humidity reference height (m) 261 rn_zu = 10. ! Wind vector reference height (m) 262 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012) 263 ln_Cd_L15 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015) 264 ln_taudif = .false. ! HF tau contribution: use "mean of stress module - module of the mean stress" data 265 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow) 266 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.) 267 rn_vfac = 0. ! multiplicative factor for ocean & ice velocity used to 268 ! ! calculate the wind stress (0.=absolute or 1.=relative winds) 269 261 rn_zqt = 10. ! Air temperature & humidity reference height (m) 262 rn_zu = 10. ! Wind vector reference height (m) 263 ln_Cd_L12 = .false. ! air-ice drags = F(ice conc.) (Lupkes et al. 2012) 264 ln_Cd_L15 = .false. ! air-ice drags = F(ice conc.) (Lupkes et al. 2015) 265 ! ! - module of the mean stress" data 266 rn_pfac = 1. ! multipl. factor for precipitation (total & snow) 267 rn_efac = 1. ! multipl. factor for evaporation (0. or 1.) 268 rn_vfac = 0. ! multipl. factor for ocean & ice velocity 269 ! ! used to calculate the wind stress 270 ! ! (0. => absolute or 1. => relative winds) 271 ln_skin_cs = .false. ! use the cool-skin parameterization 272 ln_skin_wl = .false. ! use the warm-layer parameterization 273 ! ! ==> only available in ECMWF and COARE algorithms 274 ln_humi_sph = .true. ! humidity "sn_humi" is specific humidity [kg/kg] 275 ln_humi_dpt = .false. ! humidity "sn_humi" is dew-point temperature [K] 276 ln_humi_rlh = .false. ! humidity "sn_humi" is relative humidity [%] 277 ! 270 278 cn_dir = './' ! root directory for the bulk data location 271 279 !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! … … 278 286 sn_tair = 't_10.15JUNE2009_fill' , 6. , 'T_10_MOD', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 279 287 sn_humi = 'q_10.15JUNE2009_fill' , 6. , 'Q_10_MOD', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 288 sn_hpgi = 'NONE' , 24. , 'uhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG' , '' 289 sn_hpgj = 'NONE' , 24. , 'vhpg' , .false. , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG' , '' 280 290 sn_prec = 'ncar_precip.15JUNE2009_fill', -1. , 'PRC_MOD1', .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 281 291 sn_snow = 'ncar_precip.15JUNE2009_fill', -1. , 'SNOW' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 282 292 sn_slp = 'slp.15JUNE2009_fill' , 6. , 'SLP' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 283 sn_tdif = 'taudif_core' , 24 , 'taudif' , .false. , .true. , 'yearly' , 'weights_core_orca2_bilinear_noc.nc' , '' , '' 293 / 294 !----------------------------------------------------------------------- 295 &namsbc_abl ! Atmospheric Boundary Layer formulation (ln_abl = T) 296 !----------------------------------------------------------------------- 297 cn_dir = './' ! root directory for the location of the ABL grid file 298 cn_dom = 'dom_cfg_abl.nc' 299 300 cn_ablrst_in = "restart_abl" ! suffix of abl restart name (input) 301 cn_ablrst_out = "restart_abl" ! suffix of abl restart name (output) 302 cn_ablrst_indir = "." ! directory to read input abl restarts 303 cn_ablrst_outdir = "." ! directory to write output abl restarts 304 305 ln_hpgls_frc = .false. 306 ln_geos_winds = .false. 307 nn_dyn_restore = 2 ! restoring option for dynamical ABL variables: = 0 no restoring 308 ! = 1 equatorial restoring 309 ! = 2 global restoring 310 rn_ldyn_min = 4.5 ! magnitude of the nudging on ABL dynamics at the bottom of the ABL [hour] 311 rn_ldyn_max = 1.5 ! magnitude of the nudging on ABL dynamics at the top of the ABL [hour] 312 rn_ltra_min = 4.5 ! magnitude of the nudging on ABL tracers at the bottom of the ABL [hour] 313 rn_ltra_max = 1.5 ! magnitude of the nudging on ABL tracers at the top of the ABL [hour] 314 nn_amxl = 0 ! mixing length: = 0 Deardorff 80 length-scale 315 ! = 1 length-scale based on the distance to the PBL height 316 ! = 2 Bougeault & Lacarrere 89 length-scale 317 rn_Cm = 0.0667 ! 0.126 in MesoNH 318 rn_Ct = 0.1667 ! 0.143 in MesoNH 319 rn_Ce = 0.4 ! 0.4 in MesoNH 320 rn_Ceps = 0.7 ! 0.85 in MesoNH 321 rn_Rod = 0.15 ! c0 in RMCA17 mixing length formulation (not yet implemented) 322 rn_Ric = 0.139 ! Critical Richardson number (to compute PBL height and diffusivities) 284 323 / 285 324 !----------------------------------------------------------------------- … … 375 414 nn_chldta = 0 ! RGB : Chl data (=1) or cst value (=0) 376 415 rn_si1 = 23.0 ! 2BD : longest depth of extinction 377 416 378 417 cn_dir = './' ! root directory for the chlorophyl data location 379 418 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! … … 443 482 / 444 483 !----------------------------------------------------------------------- 445 &namsbc_isf ! Top boundary layer (ISF) (ln_isfcav =T : read (ln_read_cfg=T) 484 &namsbc_isf ! Top boundary layer (ISF) (ln_isfcav =T : read (ln_read_cfg=T) 446 485 !----------------------------------------------------------------------- or set or usr_def_zgr ) 447 ! ! type of top boundary layer 486 ! ! type of top boundary layer 448 487 nn_isf = 1 ! ice shelf melting/freezing 449 ! 1 = presence of ISF ; 2 = bg03 parametrisation 488 ! 1 = presence of ISF ; 2 = bg03 parametrisation 450 489 ! 3 = rnf file for ISF ; 4 = ISF specified freshwater flux 451 490 ! options 1 and 4 need ln_isfcav = .true. (domzgr) … … 470 509 !* nn_isf = 3 case 471 510 sn_rnfisf = 'rnfisf' , -12. ,'sofwfisf' , .false. , .true. , 'yearly' , '' , '' , '' 472 !* nn_isf = 2 and 3 cases 511 !* nn_isf = 2 and 3 cases 473 512 sn_depmax_isf ='rnfisf' , -12. ,'sozisfmax', .false. , .true. , 'yearly' , '' , '' , '' 474 513 sn_depmin_isf ='rnfisf' , -12. ,'sozisfmin', .false. , .true. , 'yearly' , '' , '' , '' … … 477 516 / 478 517 !----------------------------------------------------------------------- 479 &namsbc_iscpl ! land ice / ocean coupling option (ln_isfcav =T : read (ln_read_cfg=T) 518 &namsbc_iscpl ! land ice / ocean coupling option (ln_isfcav =T : read (ln_read_cfg=T) 480 519 !----------------------------------------------------------------------- or set or usr_def_zgr ) 481 520 nn_drown = 10 ! number of iteration of the extrapolation loop (fill the new wet cells) … … 572 611 !----------------------------------------------------------------------- 573 612 ln_tide = .false. ! Activate tides 574 ln_tide_pot = . true.! use tidal potential forcing613 ln_tide_pot = .false. ! use tidal potential forcing 575 614 ln_scal_load = .false. ! Use scalar approximation for 576 615 rn_scal_load = 0.094 ! load potential 577 616 ln_read_load = .false. ! Or read load potential from file 578 617 cn_tide_load = 'tide_LOAD_grid_T.nc' ! filename for load potential 579 ! 618 ! 580 619 ln_tide_ramp = .false. ! Use linear ramp for tides at startup 581 620 rdttideramp = 0. ! ramp duration in days … … 656 695 filtide = 'bdydta/amm12_bdytide_' ! file name root of tidal forcing files 657 696 ln_bdytide_2ddta = .false. ! 658 ln_bdytide_conj = .false. ! 697 ln_bdytide_conj = .false. ! 659 698 / 660 699 … … 683 722 !----------------------------------------------------------------------- 684 723 rn_Cd0 = 1.e-3 ! drag coefficient [-] 685 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0) 724 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0) 686 725 rn_Cdmax = 0.1 ! drag value maximum [-] (logarithmic drag) 687 726 rn_ke0 = 2.5e-3 ! background kinetic energy [m2/s2] (non-linear cases) … … 694 733 !----------------------------------------------------------------------- 695 734 rn_Cd0 = 1.e-3 ! drag coefficient [-] 696 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0) 735 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0) 697 736 rn_Cdmax = 0.1 ! drag value maximum [-] (logarithmic drag) 698 737 rn_ke0 = 2.5e-3 ! background kinetic energy [m2/s2] (non-linear cases) … … 761 800 nn_cen_v = 4 ! =2/4, vertical 2nd order CEN / 4th order COMPACT 762 801 ln_traadv_fct = .false. ! FCT scheme 763 nn_fct_h = 2 ! =2/4, horizontal 2nd / 4th order 764 nn_fct_v = 2 ! =2/4, vertical 2nd / COMPACT 4th order 802 nn_fct_h = 2 ! =2/4, horizontal 2nd / 4th order 803 nn_fct_v = 2 ! =2/4, vertical 2nd / COMPACT 4th order 765 804 ln_traadv_mus = .false. ! MUSCL scheme 766 805 ln_mus_ups = .false. ! use upstream scheme near river mouths … … 783 822 ln_traldf_triad = .false. ! iso-neutral (triad operator) 784 823 ! 785 ! ! iso-neutral options: 824 ! ! iso-neutral options: 786 825 ln_traldf_msc = .false. ! Method of Stabilizing Correction (both operators) 787 826 rn_slpmax = 0.01 ! slope limit (both operators) … … 793 832 nn_aht_ijk_t = 0 ! space/time variation of eddy coefficient: 794 833 ! ! =-20 (=-30) read in eddy_diffusivity_2D.nc (..._3D.nc) file 795 ! ! = 0 constant 796 ! ! = 10 F(k) =ldf_c1d 797 ! ! = 20 F(i,j) =ldf_c2d 834 ! ! = 0 constant 835 ! ! = 10 F(k) =ldf_c1d 836 ! ! = 20 F(i,j) =ldf_c2d 798 837 ! ! = 21 F(i,j,t) =Treguier et al. JPO 1997 formulation 799 838 ! ! = 30 F(i,j,k) =ldf_c2d * ldf_c1d 800 839 ! ! = 31 F(i,j,k,t)=F(local velocity and grid-spacing) 801 ! ! time invariant coefficients: aht0 = 1/2 Ud*Ld (lap case) 840 ! ! time invariant coefficients: aht0 = 1/2 Ud*Ld (lap case) 802 841 ! ! or = 1/12 Ud*Ld^3 (blp case) 803 842 rn_Ud = 0.01 ! lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30) … … 825 864 nn_aei_ijk_t = 0 ! space/time variation of eddy coefficient: 826 865 ! ! =-20 (=-30) read in eddy_induced_velocity_2D.nc (..._3D.nc) file 827 ! ! = 0 constant 828 ! ! = 10 F(k) =ldf_c1d 829 ! ! = 20 F(i,j) =ldf_c2d 866 ! ! = 0 constant 867 ! ! = 10 F(k) =ldf_c1d 868 ! ! = 20 F(i,j) =ldf_c2d 830 869 ! ! = 21 F(i,j,t) =Treguier et al. JPO 1997 formulation 831 870 ! ! = 30 F(i,j,k) =ldf_c2d * ldf_c1d 832 ! ! time invariant coefficients: aei0 = 1/2 Ue*Le 871 ! ! time invariant coefficients: aei0 = 1/2 Ue*Le 833 872 rn_Ue = 0.02 ! lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30) 834 873 rn_Le = 200.e+3 ! lateral diffusive length [m] (nn_aht_ijk_t= 0, 10) … … 870 909 rn_lf_cutoff = 5.0 ! cutoff frequency for low-pass filter [days] 871 910 rn_zdef_max = 0.9 ! maximum fractional e3t deformation 872 ln_vvl_dbg = . true.! debug prints (T/F)911 ln_vvl_dbg = .false. ! debug prints (T/F) 873 912 / 874 913 !----------------------------------------------------------------------- … … 890 929 ln_dynvor_eeT = .false. ! energy conserving scheme (een using e3t) 891 930 ln_dynvor_een = .false. ! energy & enstrophy scheme 892 nn_een_e3f = 0 ! =0 e3f = mi(mj(e3t))/4 931 nn_een_e3f = 0 ! =0 e3f = mi(mj(e3t))/4 893 932 ! ! =1 e3f = mi(mj(e3t))/mi(mj( tmask)) 894 933 ln_dynvor_msk = .false. ! vorticity multiplied by fmask (=T) ==>>> PLEASE DO NOT ACTIVATE … … 935 974 ! ! =-30 read in eddy_viscosity_3D.nc file 936 975 ! ! =-20 read in eddy_viscosity_2D.nc file 937 ! ! = 0 constant 976 ! ! = 0 constant 938 977 ! ! = 10 F(k)=c1d 939 978 ! ! = 20 F(i,j)=F(grid spacing)=c2d … … 941 980 ! ! = 31 F(i,j,k)=F(grid spacing and local velocity) 942 981 ! ! = 32 F(i,j,k)=F(local gridscale and deformation rate) 943 ! ! time invariant coefficients : ahm = 1/2 Uv*Lv (lap case) 982 ! ! time invariant coefficients : ahm = 1/2 Uv*Lv (lap case) 944 983 ! ! or = 1/12 Uv*Lv^3 (blp case) 945 984 rn_Uv = 0.1 ! lateral viscous velocity [m/s] (nn_ahm_ijk_t= 0, 10, 20, 30) … … 1065 1104 ! = 0 constant 10 m length scale 1066 1105 ! = 1 0.5m at the equator to 30m poleward of 40 degrees 1067 rn_eice = 4 ! below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4 1106 rn_eice = 4 ! below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4 1068 1107 / 1069 1108 !----------------------------------------------------------------------- … … 1323 1362 ln_ctl = .FALSE. ! Toggle all report printing on/off (T/F); Ignored if sn_cfctl%l_config is T 1324 1363 sn_cfctl%l_config = .TRUE. ! IF .true. then control which reports are written with the following 1325 sn_cfctl%l_runstat = . FALSE.! switches and which areas produce reports with the proc integer settings.1364 sn_cfctl%l_runstat = .TRUE. ! switches and which areas produce reports with the proc integer settings. 1326 1365 sn_cfctl%l_trcstat = .FALSE. ! The default settings for the proc integers should ensure 1327 1366 sn_cfctl%l_oceout = .FALSE. ! that all areas report. -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/SPITZ12/EXPREF/namelist_cfg
r11536 r12154 107 107 sn_snow = 'MARv3.6-9km-Svalbard-2hourly_spitz' , 2. , 'snow' , .true. , .false. , 'yearly' , 'weights_bilin', '' , '' 108 108 sn_slp = 'MARv3.6-9km-Svalbard-2hourly_spitz' , 2. , 'slp' , .true. , .false. , 'yearly' , 'weights_bilin', '' , '' 109 sn_tdif = 'MARv3.6-9km-Svalbard-2hourly_spitz' , 2. , 'tdif' , .true. , .false. , 'yearly' , 'weights_bilin', '' , ''110 109 / 111 110 !----------------------------------------------------------------------- -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/ref_cfgs.txt
r9775 r12154 8 8 ORCA2_SAS_ICE OCE ICE NST SAS 9 9 ORCA2_ICE_PISCES OCE TOP ICE NST 10 ORCA2_ICE_ABL OCE ICE ABL 11 ORCA2_SAS_ICE_ABL OCE SAS ICE ABL 12 ORCA2_ICE OCE ICE 10 13 SPITZ12 OCE ICE 14 eORCA025_ICE OCE ICE 15 eORCA025_ICE_ABL OCE ICE ABL 16 eORCA025_SAS_ICE_ABL OCE SAS ICE ABL -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ABL/ablmod.F90
r11937 r12154 17 17 USE phycst ! physical constants 18 18 USE dom_oce, ONLY : tmask 19 USE sbc_oce, ONLY : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1 20 USE sbcblk ! use some physical constants for flux computation 19 USE sbc_oce, ONLY : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa 20 USE sbcblk ! use rn_?fac 21 USE sbcblk_phy ! use some physical constants for flux computation 21 22 ! 22 23 USE prtctl ! Print control (prt_ctl routine) … … 100 101 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point) 101 102 #endif 102 103 REAL(wp), DIMENSION(1:jpi,1:jpj ) :: z rhoa, zwnd_i, zwnd_j103 ! 104 REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j 104 105 REAL(wp), DIMENSION(1:jpi,2:jpka ) :: zCF 105 106 REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: z_cft !--FL--to be removed after the test phase … … 529 530 ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta ) 530 531 zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa ) 531 zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ... 532 & ( R_dry*ztemp * ( 1._wp + rctv0*zhumi ) ) 532 !zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ... 533 ! & ( R_dry*ztemp * ( 1._wp + rctv0*zhumi ) ) 534 zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 533 535 psen ( ji, jj ) = cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) 534 536 pevp ( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) ) 535 zrhoa( ji, jj ) = zcff537 rhoa( ji, jj ) = zcff 536 538 END DO 537 539 END DO … … 551 553 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 552 554 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 553 zztmp = zrhoa(ji,jj) * pcd_du(ji,jj)555 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 554 556 555 557 pwndm (ji,jj) = zcff … … 593 595 zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 594 596 595 ptaui_ice(ji,jj) = 0.5_wp * ( zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) &596 & + zrhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) &597 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 598 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 597 599 & * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 598 ptauj_ice(ji,jj) = 0.5_wp * ( zrhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) &599 & + zrhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) &600 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 601 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 600 602 & * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 601 603 END DO -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ABL/sbcabl.F90
r11858 r12154 22 22 USE sbc_oce ! Surface boundary condition: ocean fields 23 23 USE sbcblk ! Surface boundary condition: bulk formulae 24 USE sbcblk_phy ! Surface boundary condition: bulk formulae 24 25 USE dom_oce, ONLY : tmask 25 26 ! … … 93 94 IF( nn_dyn_restore < 0 .OR. nn_dyn_restore > 2 ) & 94 95 & CALL ctl_stop( 'abl_init : bad flag, nn_dyn_restore must be 0, 1 or 2 ' ) 95 ! 96 96 97 !!--------------------------------------------------------------------- 97 98 !! Control prints … … 215 216 WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 216 217 ! Check that restoring coefficients are between 0 and 1 217 !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp ) &218 !IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp ) &219 218 IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp ) & 220 219 & CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' ) 221 !IF( zcff > 1._wp .OR. zcff < 0._wp ) &222 220 IF( zcff - nn_fsbc > 0.001_wp .OR. zcff < 0._wp ) & 223 221 & CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_min' ) … … 236 234 WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 237 235 ! Check that restoring coefficients are between 0 and 1 238 !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp ) &239 236 IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp ) & 240 237 & CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' ) 241 !IF( zcff > 1._wp .OR. zcff < 0._wp ) &242 238 IF( zcff - nn_fsbc > 0.001_wp .OR. zcff < 0._wp ) & 243 239 & CALL ctl_stop( 'abl_init : wrong value for rn_ltra_min' ) … … 294 290 tke_abl(:,:,:,nt_a ) = 0._wp 295 291 ENDIF 292 293 rhoa(:,:) = rho_air( tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), sf(jp_slp)%fnow(:,:,1) ) !!GS: rhoa must be (re)computed here here to avoid division by zero in blk_ice_1 (TBI) 296 294 297 295 END SUBROUTINE sbc_abl_init … … 341 339 & tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in 342 340 & sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in 341 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in 343 342 & zssq, zcd_du, zsen, zevp ) ! =>> out 344 343 -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ICE/icesbc.F90
r11575 r12154 27 27 USE lbclnk ! lateral boundary conditions (or mpp links) 28 28 USE timing ! Timing 29 USE fldread !!GS: needed by agrif 29 30 30 31 IMPLICIT NONE … … 71 72 SELECT CASE( ksbc ) 72 73 CASE( jp_usr ) ; CALL usrdef_sbc_ice_tau( kt ) ! user defined formulation 73 CASE( jp_blk ) ; CALL blk_ice_tau ! Bulk formulation 74 CASE( jp_blk ) ; CALL blk_ice_1( sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & 75 & sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & 76 & sf(jp_slp )%fnow(:,:,1), u_ice, v_ice, tm_su , & ! inputs 77 & putaui = utau_ice, pvtaui = vtau_ice ) ! outputs 78 ! CASE( jp_abl ) utau_ice & vtau_ice are computed in ablmod 74 79 CASE( jp_purecpl ) ; CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) ! Coupled formulation 75 80 END SELECT … … 143 148 CASE( jp_usr ) !--- user defined formulation 144 149 CALL usrdef_sbc_ice_flx( kt, h_s, h_i ) 145 CASE( jp_blk ) !--- bulk formulation 146 CALL blk_ice_flx ( t_su, h_s, h_i, alb_ice ) ! 150 CASE( jp_blk, jp_abl ) !--- bulk formulation & ABL formulation 151 CALL blk_ice_2 ( t_su, h_s, h_i, alb_ice, sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & 152 & sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) ) ! 147 153 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 148 154 IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist ) -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ICE/icevar.F90
r11732 r12154 115 115 ! 116 116 ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction 117 117 ! 118 !!GS: tm_su always needed by ABL over sea-ice 119 ALLOCATE( z1_at_i(jpi,jpj) ) 120 WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:) 121 ELSEWHERE ; z1_at_i(:,:) = 0._wp 122 END WHERE 123 tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 124 WHERE( at_i(:,:)<=epsi20 ) tm_su(:,:) = rt0 125 ! 118 126 ! The following fields are calculated for diagnostics and outputs only 119 127 ! ==> Do not use them for other purposes 120 128 IF( kn > 1 ) THEN 121 129 ! 122 ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 123 WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:) 124 ELSEWHERE ; z1_at_i(:,:) = 0._wp 125 END WHERE 130 ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 126 131 WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:) 127 132 ELSEWHERE ; z1_vt_i(:,:) = 0._wp … … 136 141 ! 137 142 ! ! mean temperature (K), salinity and age 138 tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)139 143 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 140 144 om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:) … … 154 158 ! ! put rt0 where there is no ice 155 159 WHERE( at_i(:,:)<=epsi20 ) 156 tm_su(:,:) = rt0157 160 tm_si(:,:) = rt0 158 161 tm_i (:,:) = rt0 … … 165 168 END WHERE 166 169 ! 167 DEALLOCATE( z1_ at_i , z1_vt_i , z1_vt_s )170 DEALLOCATE( z1_vt_i , z1_vt_s ) 168 171 ! 169 172 ENDIF 173 ! 174 DEALLOCATE( z1_at_i ) 170 175 ! 171 176 END SUBROUTINE ice_var_agg -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DIA/diawri.F90
r11536 r12154 26 26 !!---------------------------------------------------------------------- 27 27 USE oce ! ocean dynamics and tracers 28 USE abl ! abl variables in case ln_abl = .true. 28 29 USE dom_oce ! ocean space and time domain 29 30 USE phycst ! physical constants … … 66 67 PUBLIC dia_wri_state 67 68 PUBLIC dia_wri_alloc ! Called by nemogcm module 68 69 #if ! defined key_iomput 70 PUBLIC dia_wri_alloc_abl ! Called by sbcabl module (if ln_abl = .true.) 71 #endif 69 72 INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file 70 73 INTEGER :: nb_T , ndim_bT ! grid_T file … … 72 75 INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file 73 76 INTEGER :: nid_W, nz_W, nh_W ! grid_W file 77 INTEGER :: nid_A, nz_A, nh_A, ndim_A, ndim_hA ! grid_ABL file 74 78 INTEGER :: ndex(1) ! ??? 75 79 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 80 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL 76 81 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 77 82 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT … … 414 419 & ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 415 420 ! 416 421 dia_wri_alloc = MAXVAL(ierr) 417 422 CALL mpp_sum( 'diawri', dia_wri_alloc ) 418 423 ! 419 424 END FUNCTION dia_wri_alloc 420 421 425 426 INTEGER FUNCTION dia_wri_alloc_abl() 427 !!---------------------------------------------------------------------- 428 ALLOCATE( ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl) 429 CALL mpp_sum( 'diawri', dia_wri_alloc_abl ) 430 ! 431 END FUNCTION dia_wri_alloc_abl 432 422 433 SUBROUTINE dia_wri( kt ) 423 434 !!--------------------------------------------------------------------- … … 440 451 INTEGER :: ierr ! error code return from allocation 441 452 INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers 453 INTEGER :: ipka ! ABL 442 454 INTEGER :: jn, ierror ! local integers 443 455 REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars … … 445 457 REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace 446 458 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace 459 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace 447 460 !!---------------------------------------------------------------------- 448 461 ! … … 478 491 ijmi = 1 ; ijma = jpj 479 492 ipk = jpk 493 IF(ln_abl) ipka = jpkam1 480 494 481 495 ! define time axis … … 580 594 & "m", ipk, gdepw_1d, nz_W, "down" ) 581 595 596 IF( ln_abl ) THEN 597 ! Define the ABL grid FILE ( nid_A ) 598 CALL dia_nam( clhstnam, nn_write, 'grid_ABL' ) 599 IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename 600 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit 601 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 602 & nit000-1, zjulian, rdt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set ) 603 CALL histvert( nid_A, "ght_abl", "Vertical T levels", & ! Vertical grid: gdept 604 & "m", ipka, ght_abl(2:jpka), nz_A, "up" ) 605 ! ! Index of ocean points 606 ALLOCATE( zw3d_abl(jpi,jpj,ipka) ) 607 zw3d_abl(:,:,:) = 1._wp 608 CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A ) ! volume 609 CALL wheneq( jpi*jpj , zw3d_abl, 1, 1., ndex_hA, ndim_hA ) ! surface 610 DEALLOCATE(zw3d_abl) 611 ENDIF 582 612 583 613 ! Declare all the output fields as NETCDF variables … … 629 659 CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm 630 660 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 631 ! 661 ! 662 IF( ln_abl ) THEN 663 CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl 664 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 665 CALL histdef( nid_A, "q_abl", "Humidity" , "kg/kg" , & ! q_abl 666 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 667 CALL histdef( nid_A, "u_abl", "Atmospheric U-wind " , "m/s" , & ! u_abl 668 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 669 CALL histdef( nid_A, "v_abl", "Atmospheric V-wind " , "m/s" , & ! v_abl 670 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 671 CALL histdef( nid_A, "tke_abl", "Atmospheric TKE " , "m2/s2" , & ! tke_abl 672 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 673 CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s" , & ! avm_abl 674 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 675 CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2", & ! avt_abl 676 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 677 CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height " , "m", & ! pblh 678 & jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 679 #if defined key_si3 680 CALL histdef( nid_A, "oce_frac", "Fraction of open ocean" , " ", & ! ato_i 681 & jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 682 #endif 683 CALL histend( nid_A, snc4chunks=snc4set ) 684 ENDIF 685 ! 632 686 IF( ln_icebergs ) THEN 633 687 CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , & … … 787 841 CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction 788 842 CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed 789 ! 843 ! 844 IF( ln_abl ) THEN 845 ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 846 IF( ln_mskland ) THEN 847 DO jk=1,jpka 848 zw3d_abl(:,:,jk) = tmask(:,:,1) 849 END DO 850 ELSE 851 zw3d_abl(:,:,:) = 1._wp 852 ENDIF 853 CALL histwrite( nid_A, "pblh" , it, pblh(:,:) *zw3d_abl(:,:,1 ), ndim_hA, ndex_hA ) ! pblh 854 CALL histwrite( nid_A, "u_abl" , it, u_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! u_abl 855 CALL histwrite( nid_A, "v_abl" , it, v_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! v_abl 856 CALL histwrite( nid_A, "t_abl" , it, tq_abl (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! t_abl 857 CALL histwrite( nid_A, "q_abl" , it, tq_abl (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! q_abl 858 CALL histwrite( nid_A, "tke_abl", it, tke_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! tke_abl 859 CALL histwrite( nid_A, "avm_abl", it, avm_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avm_abl 860 CALL histwrite( nid_A, "avt_abl", it, avt_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avt_abl 861 #if defined key_si3 862 CALL histwrite( nid_A, "oce_frac" , it, ato_i(:,:) , ndim_hA, ndex_hA ) ! ato_i 863 #endif 864 DEALLOCATE(zw3d_abl) 865 ENDIF 866 ! 790 867 IF( ln_icebergs ) THEN 791 868 ! … … 857 934 CALL histclo( nid_V ) 858 935 CALL histclo( nid_W ) 936 IF(ln_abl) CALL histclo( nid_A ) 859 937 ENDIF 860 938 ! … … 926 1004 CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity 927 1005 ENDIF 1006 IF ( ln_abl ) THEN 1007 CALL iom_rstput ( 0, 0, inum, "uz1_abl", u_abl(:,:,2,nt_a ) ) ! now first level i-wind 1008 CALL iom_rstput ( 0, 0, inum, "vz1_abl", v_abl(:,:,2,nt_a ) ) ! now first level j-wind 1009 CALL iom_rstput ( 0, 0, inum, "tz1_abl", tq_abl(:,:,2,nt_a,1) ) ! now first level temperature 1010 CALL iom_rstput ( 0, 0, inum, "qz1_abl", tq_abl(:,:,2,nt_a,2) ) ! now first level humidity 1011 ENDIF 928 1012 929 1013 #if defined key_si3 -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/IOM/in_out_manager.F90
r11536 r12154 87 87 LOGICAL :: lrst_oce !: logical to control the oce restart write 88 88 LOGICAL :: lrst_ice !: logical to control the ice restart write 89 LOGICAL :: lrst_abl !: logical to control the abl restart write 89 90 INTEGER :: numror = 0 !: logical unit for ocean restart (read). Init to 0 is needed for SAS (in daymod.F90) 90 91 INTEGER :: numrir !: logical unit for ice restart (read) 92 INTEGER :: numrar !: logical unit for abl restart (read) 91 93 INTEGER :: numrow !: logical unit for ocean restart (write) 92 94 INTEGER :: numriw !: logical unit for ice restart (write) 95 INTEGER :: numraw !: logical unit for abl restart (write) 93 96 INTEGER :: nrst_lst !: number of restart to output next 94 97 -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/IOM/iom.F90
r12109 r12154 29 29 USE lib_mpp ! MPP library 30 30 #if defined key_iomput 31 USE sbc_oce , ONLY : nn_fsbc ! ocean space and time domain31 USE sbc_oce , ONLY : nn_fsbc, ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1 32 32 USE trc_oce , ONLY : nn_dttrc ! !: frequency of step on passive tracers 33 33 USE icb_oce , ONLY : nclasses, class_num ! !: iceberg classes … … 113 113 ! 114 114 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zt_bnds, zw_bnds 115 REAL(wp), DIMENSION(2,jpkam1) :: za_bnds ! ABL vertical boundaries 115 116 LOGICAL :: ll_tmppatch = .TRUE. !: seb: patch before we remove periodicity 116 117 INTEGER :: nldi_save, nlei_save !: and close boundaries in output files … … 200 201 ! vertical grid definition 201 202 IF(.NOT.llrst_context) THEN 202 CALL iom_set_axis_attr( "deptht", paxis = gdept_1d ) 203 CALL iom_set_axis_attr( "depthu", paxis = gdept_1d ) 204 CALL iom_set_axis_attr( "depthv", paxis = gdept_1d ) 205 CALL iom_set_axis_attr( "depthw", paxis = gdepw_1d ) 206 203 CALL iom_set_axis_attr( "deptht", paxis = gdept_1d ) 204 CALL iom_set_axis_attr( "depthu", paxis = gdept_1d ) 205 CALL iom_set_axis_attr( "depthv", paxis = gdept_1d ) 206 CALL iom_set_axis_attr( "depthw", paxis = gdepw_1d ) 207 208 ! ABL 209 IF( .NOT. ALLOCATED(ght_abl) ) THEN ! force definition for xml files (xios) 210 ALLOCATE( ght_abl(jpka), ghw_abl(jpka), e3t_abl(jpka), e3w_abl(jpka) ) ! default allocation needed by iom 211 ght_abl(:) = -1._wp ; ghw_abl(:) = -1._wp 212 e3t_abl(:) = -1._wp ; e3w_abl(:) = -1._wp 213 ENDIF 214 CALL iom_set_axis_attr( "ght_abl", ght_abl(2:jpka) ) 215 CALL iom_set_axis_attr( "ghw_abl", ghw_abl(2:jpka) ) 216 207 217 ! Add vertical grid bounds 208 218 jkmin = MIN(2,jpk) ! in case jpk=1 (i.e. sas2D) … … 213 223 zw_bnds(2,1:jpkm1 ) = gdepw_1d(jkmin:jpk) 214 224 zw_bnds(2,jpk: ) = gdepw_1d(jpk) + e3t_1d(jpk) 215 CALL iom_set_axis_attr( "deptht", bounds=zw_bnds ) 216 CALL iom_set_axis_attr( "depthu", bounds=zw_bnds ) 217 CALL iom_set_axis_attr( "depthv", bounds=zw_bnds ) 218 CALL iom_set_axis_attr( "depthw", bounds=zt_bnds ) 225 CALL iom_set_axis_attr( "deptht", bounds=zw_bnds ) 226 CALL iom_set_axis_attr( "depthu", bounds=zw_bnds ) 227 CALL iom_set_axis_attr( "depthv", bounds=zw_bnds ) 228 CALL iom_set_axis_attr( "depthw", bounds=zt_bnds ) 229 230 ! ABL 231 za_bnds(1,:) = ghw_abl(1:jpkam1) 232 za_bnds(2,:) = ghw_abl(2:jpka ) 233 CALL iom_set_axis_attr( "ght_abl", bounds=za_bnds ) 234 za_bnds(1,:) = ght_abl(2:jpka ) 235 za_bnds(2,:) = ght_abl(2:jpka ) + e3w_abl(2:jpka) 236 CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds ) 237 219 238 CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) ) 220 239 # if defined key_si3 … … 1147 1166 WRITE(cldmspc , fmt='(i1)') idmspc 1148 1167 ! 1149 IF( idmspc < irankpv ) THEN 1150 CALL ctl_stop( TRIM(clinfo), 'The file has only '//cldmspc//' spatial dimension', & 1151 & 'it is impossible to read a '//clrankpv//'D array from this file...' ) 1152 ELSEIF( idmspc == irankpv ) THEN 1168 !!GS: we consider 2D data as 3D data with vertical dim size = 1 1169 !IF( idmspc < irankpv ) THEN 1170 ! CALL ctl_stop( TRIM(clinfo), 'The file has only '//cldmspc//' spatial dimension', & 1171 ! & 'it is impossible to read a '//clrankpv//'D array from this file...' ) 1172 !ELSEIF( idmspc == irankpv ) THEN 1173 IF( idmspc == irankpv ) THEN 1153 1174 IF( PRESENT(pv_r1d) .AND. idom /= jpdom_unknown ) & 1154 1175 & CALL ctl_stop( TRIM(clinfo), 'case not coded...You must use jpdom_unknown' ) … … 1972 1993 ! 1973 1994 INTEGER :: ji, jj, jn, ni, nj 1974 INTEGER :: icnr, jcnr 1975 ! 1995 INTEGER :: icnr, jcnr ! Offset such that the vertex coordinate (i+icnr,j+jcnr) 1996 ! ! represents the bottom-left corner of cell (i,j) 1976 1997 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z_bnds ! Lat/lon coordinates of the vertices of cell (i,j) 1977 1998 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_fld ! Working array to determine where to rotate cells … … 2143 2164 f_op%timestep = nn_fsbc ; f_of%timestep = 0 ; CALL iom_set_field_attr('SBC' , freq_op=f_op, freq_offset=f_of) 2144 2165 f_op%timestep = nn_fsbc ; f_of%timestep = 0 ; CALL iom_set_field_attr('SBC_scalar' , freq_op=f_op, freq_offset=f_of) 2166 f_op%timestep = nn_fsbc ; f_of%timestep = 0 ; CALL iom_set_field_attr('ABL' , freq_op=f_op, freq_offset=f_of) 2145 2167 f_op%timestep = nn_dttrc ; f_of%timestep = 0 ; CALL iom_set_field_attr('ptrc_T' , freq_op=f_op, freq_offset=f_of) 2146 2168 f_op%timestep = nn_dttrc ; f_of%timestep = 0 ; CALL iom_set_field_attr('diad_T' , freq_op=f_op, freq_offset=f_of) -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/IOM/iom_nf90.F90
r11536 r12154 19 19 !!---------------------------------------------------------------------- 20 20 USE dom_oce ! ocean space and time domain 21 USE sbc_oce, ONLY: jpka, ght_abl ! abl vertical level number and height 21 22 USE lbclnk ! lateal boundary condition / mpp exchanges 22 23 USE iom_def ! iom variables definitions … … 56 57 LOGICAL , INTENT(in ) :: ldok ! check the existence 57 58 INTEGER, DIMENSION(2,5), INTENT(in ), OPTIONAL :: kdompar ! domain parameters: 58 INTEGER , INTENT(in ), OPTIONAL :: kdlev ! size of the third dimension59 INTEGER , INTENT(in ), OPTIONAL :: kdlev ! size of the ice/abl third dimension 59 60 60 61 CHARACTER(LEN=256) :: clinfo ! info character … … 69 70 INTEGER :: ihdf5 ! local variable for retrieval of value for NF90_HDF5 70 71 LOGICAL :: llclobber ! local definition of ln_clobber 71 INTEGER :: ilevels 72 INTEGER :: ilevels ! vertical levels 72 73 !--------------------------------------------------------------------- 73 74 ! … … 76 77 ! 77 78 ! !number of vertical levels 78 IF( PRESENT(kdlev) ) THEN ; ilevels = kdlev ! use input value (useful for sea-ice)79 ELSE ; ilevels = jpk ! by default jpk79 IF( PRESENT(kdlev) ) THEN ; ilevels = kdlev ! use input value (useful for sea-ice and abl) 80 ELSE ; ilevels = jpk ! by default jpk 80 81 ENDIF 81 82 ! … … 126 127 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'x', kdompar(1,1), idmy ), clinfo) 127 128 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'y', kdompar(2,1), idmy ), clinfo) 128 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', jpk, idmy ), clinfo) 129 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo) 130 IF( PRESENT(kdlev) ) & 131 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'numcat', kdlev, idmy ), clinfo) 129 IF( PRESENT(kdlev) ) THEN 130 IF( kdlev == jpka ) THEN 131 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', kdlev, idmy ), clinfo) 132 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo) 133 ELSE 134 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', jpk, idmy ), clinfo) 135 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo) 136 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'numcat', kdlev, idmy ), clinfo) 137 ENDIF 138 ELSE 139 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'nav_lev', jpk, idmy ), clinfo) 140 CALL iom_nf90_check(NF90_DEF_DIM( if90id, 'time_counter', NF90_UNLIMITED, idmy ), clinfo) 141 ENDIF 132 142 ! global attributes 133 143 CALL iom_nf90_check(NF90_PUT_ATT( if90id, NF90_GLOBAL, 'DOMAIN_number_total' , jpnij ), clinfo) … … 196 206 CHARACTER(len=*) , INTENT(in ) :: cdvar ! name of the variable 197 207 INTEGER , INTENT(in ) :: kiv ! 198 INTEGER, DIMENSION(:), INTENT( out), OPTIONAL :: kdimsz ! size of the dimensions199 INTEGER , INTENT( out), OPTIONAL :: kndims ! size of thedimensions208 INTEGER, DIMENSION(:), INTENT( out), OPTIONAL :: kdimsz ! size of each dimension 209 INTEGER , INTENT( out), OPTIONAL :: kndims ! number of dimensions 200 210 LOGICAL , INTENT( out), OPTIONAL :: lduld ! true if the last dimension is unlimited (time) 201 211 ! … … 584 594 IF( PRESENT(pv_r0d) ) THEN ; idims = 0 585 595 ELSEIF( PRESENT(pv_r1d) ) THEN 586 IF( SIZE(pv_r1d,1) == jpk) THEN ; idim3 = 3587 ELSE ; idim3 = 5596 IF(( SIZE(pv_r1d,1) == jpk ).OR.( SIZE(pv_r1d,1) == jpka )) THEN ; idim3 = 3 597 ELSE ; idim3 = 5 588 598 ENDIF 589 599 idims = 2 ; idimid(1:idims) = (/idim3,4/) 590 600 ELSEIF( PRESENT(pv_r2d) ) THEN ; idims = 3 ; idimid(1:idims) = (/1,2 ,4/) 591 601 ELSEIF( PRESENT(pv_r3d) ) THEN 592 IF( SIZE(pv_r3d,3) == jpk) THEN ; idim3 = 3593 ELSE ; idim3 = 5602 IF(( SIZE(pv_r3d,3) == jpk ).OR.( SIZE(pv_r3d,3) == jpka )) THEN ; idim3 = 3 603 ELSE ; idim3 = 5 594 604 ENDIF 595 605 idims = 4 ; idimid(1:idims) = (/1,2,idim3,4/) … … 674 684 CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, gphit(ix1:ix2, iy1:iy2) ), clinfo ) 675 685 CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lev' , idmy ), clinfo ) 676 CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, gdept_1d ), clinfo ) 686 IF (iom_file(kiomid)%nlev == jpka) THEN ; CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, ght_abl), clinfo ) 687 ELSE ; CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, gdept_1d), clinfo ) 688 ENDIF 677 689 IF( NF90_INQ_VARID( if90id, 'numcat', idmy ) == nf90_noerr ) THEN 678 690 CALL iom_nf90_check( NF90_PUT_VAR ( if90id, idmy, (/ (idlv, idlv = 1,iom_file(kiomid)%nlev) /)), clinfo ) -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/cpl_oasis3.F90
r10582 r12154 114 114 !------------------------------------------------------------------ 115 115 CALL oasis_init_comp ( ncomp_id, TRIM(cd_modname), nerror ) 116 IF 116 IF( nerror /= OASIS_Ok ) & 117 117 CALL oasis_abort (ncomp_id, 'cpl_init', 'Failure in oasis_init_comp') 118 118 … … 122 122 123 123 CALL oasis_get_localcomm ( kl_comm, nerror ) 124 IF 124 IF( nerror /= OASIS_Ok ) & 125 125 CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' ) 126 126 ! … … 149 149 150 150 ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 151 IF 151 IF( ltmp_wapatch ) THEN 152 152 nldi_save = nldi ; nlei_save = nlei 153 153 nldj_save = nldj ; nlej_save = nlej … … 217 217 ! 218 218 DO ji = 1, ksnd 219 IF 219 IF( ssnd(ji)%laction ) THEN 220 220 221 221 IF( ssnd(ji)%nct > nmaxcat ) THEN … … 228 228 DO jm = 1, kcplmodel 229 229 230 IF 230 IF( ssnd(ji)%nct .GT. 1 ) THEN 231 231 WRITE(cli2,'(i2.2)') jc 232 232 zclname = TRIM(ssnd(ji)%clname)//'_cat'//cli2 … … 234 234 zclname = ssnd(ji)%clname 235 235 ENDIF 236 IF 236 IF( kcplmodel > 1 ) THEN 237 237 WRITE(cli2,'(i2.2)') jm 238 238 zclname = 'model'//cli2//'_'//TRIM(zclname) … … 241 241 IF( agrif_fixed() /= 0 ) THEN 242 242 zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname) 243 END 243 ENDIF 244 244 #endif 245 245 IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out 246 246 CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part , (/ 2, 1 /), & 247 247 & OASIS_Out , ishape , OASIS_REAL, nerror ) 248 IF 248 IF( nerror /= OASIS_Ok ) THEN 249 249 WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname) 250 250 CALL oasis_abort ( ssnd(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' ) … … 262 262 ! 263 263 DO ji = 1, krcv 264 IF 264 IF( srcv(ji)%laction ) THEN 265 265 266 266 IF( srcv(ji)%nct > nmaxcat ) THEN … … 273 273 DO jm = 1, kcplmodel 274 274 275 IF 275 IF( srcv(ji)%nct .GT. 1 ) THEN 276 276 WRITE(cli2,'(i2.2)') jc 277 277 zclname = TRIM(srcv(ji)%clname)//'_cat'//cli2 … … 279 279 zclname = srcv(ji)%clname 280 280 ENDIF 281 IF 281 IF( kcplmodel > 1 ) THEN 282 282 WRITE(cli2,'(i2.2)') jm 283 283 zclname = 'model'//cli2//'_'//TRIM(zclname) … … 286 286 IF( agrif_fixed() /= 0 ) THEN 287 287 zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname) 288 END 288 ENDIF 289 289 #endif 290 290 IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In 291 291 CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part , (/ 2, 1 /), & 292 292 & OASIS_In , ishape , OASIS_REAL, nerror ) 293 IF 293 IF( nerror /= OASIS_Ok ) THEN 294 294 WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname) 295 295 CALL oasis_abort ( srcv(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' ) … … 310 310 IF( nerror /= OASIS_Ok ) CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef') 311 311 ! 312 IF 312 IF( ltmp_wapatch ) THEN 313 313 nldi = nldi_save ; nlei = nlei_save 314 314 nldj = nldj_save ; nlej = nlej_save … … 332 332 !!-------------------------------------------------------------------- 333 333 ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 334 IF 334 IF( ltmp_wapatch ) THEN 335 335 nldi_save = nldi ; nlei_save = nlei 336 336 nldj_save = nldj ; nlej_save = nlej … … 349 349 CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo ) 350 350 351 IF 352 IF 351 IF( ln_ctl ) THEN 352 IF( kinfo == OASIS_Sent .OR. kinfo == OASIS_ToRest .OR. & 353 353 & kinfo == OASIS_SentOut .OR. kinfo == OASIS_ToRestOut ) THEN 354 354 WRITE(numout,*) '****************' … … 368 368 ENDDO 369 369 ENDDO 370 IF 370 IF( ltmp_wapatch ) THEN 371 371 nldi = nldi_save ; nlei = nlei_save 372 372 nldj = nldj_save ; nlej = nlej_save … … 393 393 !!-------------------------------------------------------------------- 394 394 ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 395 IF 395 IF( ltmp_wapatch ) THEN 396 396 nldi_save = nldi ; nlei_save = nlei 397 397 nldj_save = nldj ; nlej_save = nlej … … 403 403 ! 404 404 DO jc = 1, srcv(kid)%nct 405 IF 405 IF( ltmp_wapatch ) THEN 406 406 IF( nimpp == 1 ) nldi = 1 407 407 IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi … … 420 420 & kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut 421 421 422 IF 422 IF( ln_ctl ) WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm) 423 423 424 IF 424 IF( llaction ) THEN 425 425 426 426 kinfo = OASIS_Rcv … … 432 432 ENDIF 433 433 434 IF 434 IF( ln_ctl ) THEN 435 435 WRITE(numout,*) '****************' 436 436 WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname … … 450 450 ENDDO 451 451 452 IF 452 IF( ltmp_wapatch ) THEN 453 453 nldi = nldi_save ; nlei = nlei_save 454 454 nldj = nldj_save ; nlej = nlej_save … … 483 483 ! 484 484 DO ji = 1, nsnd 485 IF 485 IF(ssnd(ji)%laction ) THEN 486 486 DO jm = 1, ncplmodel 487 487 IF( ssnd(ji)%nid(1,jm) /= -1 ) THEN … … 495 495 ENDDO 496 496 DO ji = 1, nrcv 497 IF 497 IF(srcv(ji)%laction ) THEN 498 498 DO jm = 1, ncplmodel 499 499 IF( srcv(ji)%nid(1,jm) /= -1 ) THEN … … 529 529 ! 530 530 DEALLOCATE( exfld ) 531 IF 531 IF(nstop == 0) THEN 532 532 CALL oasis_terminate( nerror ) 533 533 ELSE -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/cyclone.F90
r10068 r12154 137 137 zhemi = SIGN( 1. , zrlat ) 138 138 zinfl = 15.* rad ! clim inflow angle in Tropical Cyclones 139 IF 139 IF( vortex == 0 ) THEN 140 140 141 141 ! Vortex Holland reconstruct wind at each lon-lat position … … 157 157 & + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 158 158 159 IF 159 IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius 160 160 ! shape of the wind profile 161 161 zztmp = ( zrmw / ( zdist + 1.e-12 ) )**zb 162 162 zztmp = zvmax * SQRT( zztmp * EXP(1. - zztmp) ) 163 163 164 IF 164 IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 165 165 zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 166 166 ENDIF 167 167 168 168 ! !!! KILL EQ WINDS 169 ! IF 169 ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN 170 170 ! zztmp = 0. ! winds in other hemisphere 171 ! IF 172 ! ENDIF 173 ! IF 171 ! IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0. ! kill between 5N-5S 172 ! ENDIF 173 ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 174 174 ! zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) ) 175 175 ! !linear to zero between 10 and 5 … … 177 177 ! !!! / KILL EQ 178 178 179 IF 179 IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 180 180 181 181 zwnd_t = COS( zinfl ) * zztmp … … 196 196 END DO 197 197 198 ELSE IF 198 ELSE IF( vortex == 1 ) THEN 199 199 200 200 ! Vortex Willoughby reconstruct wind at each lon-lat position … … 206 206 zn = 2.1340 + 0.0077*zvmax - 0.4522*LOG(zrmw/1000.) - 0.0038*ABS( ztct(jtc,jp_lat) ) 207 207 zA = 0.5913 + 0.0029*zvmax - 0.1361*LOG(zrmw/1000.) - 0.0042*ABS( ztct(jtc,jp_lat) ) 208 IF 208 IF(zA < 0) THEN 209 209 zA=0 210 210 ENDIF … … 218 218 & + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 219 219 220 IF 220 IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius 221 221 222 222 ! shape of the wind profile 223 IF 223 IF(zdist <= zrmw) THEN ! inside the Radius of Maximum Wind 224 224 zztmp = zvmax * (zdist/zrmw)**zn 225 225 ELSE … … 227 227 ENDIF 228 228 229 IF 229 IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 230 230 zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 231 231 ENDIF 232 232 233 233 ! !!! KILL EQ WINDS 234 ! IF 234 ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN 235 235 ! zztmp = 0. ! winds in other hemisphere 236 ! IF 237 ! ENDIF 238 ! IF 236 ! IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0. ! kill between 5N-5S 237 ! ENDIF 238 ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 239 239 ! zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) ) 240 240 ! !linear to zero between 10 and 5 … … 242 242 ! !!! / KILL EQ 243 243 244 IF 244 IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 245 245 246 246 zwnd_t = COS( zinfl ) * zztmp -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/fldread.F90
r11536 r12154 167 167 IF( PRESENT(kit) ) ll_firstcall = ll_firstcall .and. kit == 1 168 168 169 IF 169 IF( nn_components == jp_iam_sas ) THEN ; it_offset = nn_fsbc 170 170 ELSE ; it_offset = 0 171 171 ENDIF … … 389 389 ENDIF 390 390 ! 391 IF 391 IF( sdjf%cltype(1:4) == 'week' ) THEN 392 392 isec_week = isec_week + ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week 393 393 llprevmth = isec_week > nsec_month ! longer time since the beginning of the week than the month … … 464 464 ENDIF 465 465 ! 466 IF 466 IF( nn_components == jp_iam_sas ) THEN ; it_offset = nn_fsbc 467 467 ELSE ; it_offset = 0 468 468 ENDIF … … 656 656 ENDIF 657 657 CASE DEFAULT 658 IF 658 IF(lk_c1d .AND. lmoor ) THEN 659 659 IF( sdjf%ln_tint ) THEN 660 660 CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) ) … … 1071 1071 imonth = kmonth 1072 1072 iday = kday 1073 IF 1073 IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week 1074 1074 isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 ) 1075 1075 llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month … … 1080 1080 ENDIF 1081 1081 ELSE ! use current day values 1082 IF 1082 IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week 1083 1083 isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week 1084 1084 llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month … … 1318 1318 1319 1319 !! get dimensions 1320 IF ( SIZE(sd%fnow, 3) > 1 ) THEN 1320 !!GS: we consider 2D data as 3D data with vertical dim size = 1 1321 !IF( SIZE(sd%fnow, 3) > 1 ) THEN 1322 IF( SIZE(sd%fnow, 3) > 0 ) THEN 1321 1323 ALLOCATE( ddims(4) ) 1322 1324 ELSE … … 1331 1333 1332 1334 CALL iom_open ( sd%wgtname, inum ) ! interpolation weights 1333 IF 1335 IF( inum > 0 ) THEN 1334 1336 1335 1337 !! determine whether we have an east-west cyclic grid … … 1640 1642 1641 1643 ref_wgts(kw)%fly_dta(:,:,:) = 0.0 1642 SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) ) 1643 CASE(1) 1644 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn) 1645 CASE DEFAULT 1644 !!GS: we consider 2D data as 3D data with vertical dim size = 1 1645 !SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) ) 1646 !CASE(1) 1647 ! CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn) 1648 !CASE DEFAULT 1646 1649 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn) 1647 END SELECT1650 !END SELECT 1648 1651 ENDIF 1649 1652 … … 1663 1666 END DO 1664 1667 1665 IF 1668 IF(ref_wgts(kw)%numwgt .EQ. 16) THEN 1666 1669 1667 1670 !! fix up halo points that we couldnt read from file … … 1689 1692 IF( jpi1 == 2 ) THEN 1690 1693 rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap 1691 SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 1692 CASE(1) 1693 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 1694 CASE DEFAULT 1694 !!GS: we consider 2D data as 3D data with vertical dim size = 1 1695 !SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 1696 !CASE(1) 1697 ! CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 1698 !CASE DEFAULT 1695 1699 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 1696 END SELECT1700 !END SELECT 1697 1701 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1698 1702 ENDIF 1699 1703 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 1700 1704 rec1(1) = 1 + ref_wgts(kw)%overlap 1701 SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 1702 CASE(1) 1703 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 1704 CASE DEFAULT 1705 !!GS: we consider 2D data as 3D data with vertical dim size = 1 1706 !SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 1707 !CASE(1) 1708 ! CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 1709 !CASE DEFAULT 1705 1710 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 1706 END SELECT1711 !END SELECT 1707 1712 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1708 1713 ENDIF … … 1746 1751 END DO 1747 1752 ! 1748 END 1753 ENDIF 1749 1754 ! 1750 1755 END SUBROUTINE fld_interp -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbc_oce.F90
r10882 r12154 11 11 !! 4.0 ! 2012-05 (C. Rousset) add attenuation coef for use in ice model 12 12 !! 4.0 ! 2016-06 (L. Brodeau) new unified bulk routine (based on AeroBulk) 13 !! 4.0 ! 2019-03 (F. Lemarié, G. Samson) add compatibility with ABL mode 13 14 !!---------------------------------------------------------------------- 14 15 … … 34 35 LOGICAL , PUBLIC :: ln_flx !: flux formulation 35 36 LOGICAL , PUBLIC :: ln_blk !: bulk formulation 37 LOGICAL , PUBLIC :: ln_abl !: Atmospheric boundary layer model 36 38 #if defined key_oasis3 37 39 LOGICAL , PUBLIC :: lk_oasis = .TRUE. !: OASIS used … … 77 79 INTEGER , PUBLIC, PARAMETER :: jp_flx = 2 !: flux formulation 78 80 INTEGER , PUBLIC, PARAMETER :: jp_blk = 3 !: bulk formulation 79 INTEGER , PUBLIC, PARAMETER :: jp_purecpl = 4 !: Pure ocean-atmosphere Coupled formulation 80 INTEGER , PUBLIC, PARAMETER :: jp_none = 5 !: for OPA when doing coupling via SAS module 81 INTEGER , PUBLIC, PARAMETER :: jp_abl = 4 !: Atmospheric boundary layer formulation 82 INTEGER , PUBLIC, PARAMETER :: jp_purecpl = 5 !: Pure ocean-atmosphere Coupled formulation 83 INTEGER , PUBLIC, PARAMETER :: jp_none = 6 !: for OPA when doing coupling via SAS module 81 84 82 85 !!---------------------------------------------------------------------- … … 107 110 INTEGER , PUBLIC :: ncpl_qsr_freq !: qsr coupling frequency per days from atmosphere 108 111 ! 109 LOGICAL , PUBLIC :: lhftau = .FALSE. !: HF tau used in TKE: mean(stress module) - module(mean stress)110 112 !! !! now ! before !! 111 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau , utau_b !: sea surface i-stress (ocean referential) [N/m2] 112 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau , vtau_b !: sea surface j-stress (ocean referential) [N/m2] 113 115 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: taum !: module of sea surface stress (at T-point) [N/m2] 114 !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads116 !! wndm is used compute surface gases exchanges in ice-free ocean or leads 115 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s] 118 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhoa !: air density at "rn_zu" m above the sea [kg/m3] !LB 116 119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr !: sea heat flux: solar [W/m2] 117 120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns , qns_b !: sea heat flux: non solar [W/m2] … … 137 140 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 138 141 142 !!--------------------------------------------------------------------- 143 !! ABL Vertical Domain size 144 !!--------------------------------------------------------------------- 145 INTEGER , PUBLIC :: jpka = 2 !: ABL number of vertical levels (default definition) 146 INTEGER , PUBLIC :: jpkam1 = 1 !: jpka-1 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ght_abl, ghw_abl !: ABL geopotential height (needed for iom) 148 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_abl, e3w_abl !: ABL vertical scale factors (needed for iom) 149 139 150 !!---------------------------------------------------------------------- 140 151 !! Sea Surface Mean fields … … 167 178 ! 168 179 ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) , & 169 & vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , STAT=ierr(1) )180 & vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) ) 170 181 ! 171 182 ALLOCATE( qns_tot(jpi,jpj) , qns (jpi,jpj) , qns_b(jpi,jpj), & … … 182 193 & ssu_m (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) , & 183 194 & ssv_m (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 184 195 ! 185 196 ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 186 197 ! 187 198 sbc_oce_alloc = MAXVAL( ierr ) 188 199 CALL mpp_sum ( 'sbc_oce', sbc_oce_alloc ) -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcapr.F90
r11536 r12154 103 103 ! 104 104 ! !* control check 105 IF 105 IF( ln_apr_obc ) THEN 106 106 IF(lwp) WRITE(numout,*) ' Inverse barometer added to OBC ssh data' 107 107 ENDIF -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk.F90
r12109 r12154 15 15 !! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk 16 16 !! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 17 !! ! ==> based on AeroBulk (http ://aerobulk.sourceforge.net/)17 !! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/) 18 18 !! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 20 !! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE) 20 21 !!---------------------------------------------------------------------- 21 22 … … 23 24 !! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition 24 25 !! sbc_blk : bulk formulation as ocean surface boundary condition 25 !! blk_oce : computes momentum, heat and freshwater fluxes over ocean 26 !! rho_air : density of (moist) air (depends on T_air, q_air and SLP 27 !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) 28 !! q_sat : saturation humidity as a function of SLP and temperature 29 !! L_vap : latent heat of vaporization of water as a function of temperature 30 !! sea-ice case only : 31 !! blk_ice_tau : provide the air-ice stress 32 !! blk_ice_flx : provide the heat and mass fluxes at air-ice interface 26 !! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE) 27 !! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE) 28 !! sea-ice case only : 29 !! blk_ice_1 : provide the air-ice stress 30 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 33 31 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 34 32 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 35 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 33 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 36 34 !!---------------------------------------------------------------------- 37 35 USE oce ! ocean dynamics and tracers … … 46 44 USE lib_fortran ! to use key_nosignedzero 47 45 #if defined key_si3 48 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif46 USE ice , ONLY : jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif 49 47 USE icethd_dh ! for CALL ice_thd_snwblow 50 48 #endif 51 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 52 USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003)53 USE sbcblk_algo_coare3p 5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)54 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31)49 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 50 USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 51 USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 52 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 45r1) 55 53 ! 56 54 USE iom ! I/O manager library … … 60 58 USE prtctl ! Print control 61 59 60 USE sbcblk_phy ! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc... 61 62 62 63 IMPLICIT NONE 63 64 PRIVATE … … 65 66 PUBLIC sbc_blk_init ! called in sbcmod 66 67 PUBLIC sbc_blk ! called in sbcmod 68 PUBLIC blk_oce_1 ! called in sbcabl 69 PUBLIC blk_oce_2 ! called in sbcabl 67 70 #if defined key_si3 68 PUBLIC blk_ice_ tau! routine called in icesbc69 PUBLIC blk_ice_ flx! routine called in icesbc71 PUBLIC blk_ice_1 ! routine called in icesbc 72 PUBLIC blk_ice_2 ! routine called in icesbc 70 73 PUBLIC blk_ice_qcn ! routine called in icesbc 71 #endif 72 73 !!Lolo: should ultimately be moved in the module with all physical constants ? 74 !!gm : In principle, yes. 75 REAL(wp), PARAMETER :: Cp_dry = 1005.0 !: Specic heat of dry air, constant pressure [J/K/kg] 76 REAL(wp), PARAMETER :: Cp_vap = 1860.0 !: Specic heat of water vapor, constant pressure [J/K/kg] 77 REAL(wp), PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg] 78 REAL(wp), PARAMETER :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg] 79 REAL(wp), PARAMETER :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 80 REAL(wp), PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 81 82 INTEGER , PARAMETER :: jpfld =10 ! maximum number of files to read 83 INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point 84 INTEGER , PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point 85 INTEGER , PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin) 86 INTEGER , PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 87 INTEGER , PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2) 88 INTEGER , PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2) 89 INTEGER , PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 90 INTEGER , PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 91 INTEGER , PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa) 92 INTEGER , PARAMETER :: jp_tdif =10 ! index of tau diff associated to HF tau (N/m2) at T-point 93 94 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read) 95 96 ! !!! Bulk parameters 97 REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...) 98 REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation 99 REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant 100 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice 101 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 102 ! 74 #endif 75 76 INTEGER , PUBLIC :: jpfld ! maximum number of files to read 77 INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point 78 INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point 79 INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin) 80 INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 81 INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2) 82 INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2) 83 INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 84 INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 85 INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa) 86 INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 87 INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 88 89 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read) 90 103 91 ! !!* Namelist namsbc_blk : bulk parameters 104 92 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) 105 93 LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003) 106 LOGICAL :: ln_COARE_3p 5 ! "COARE 3.5" algorithm (Edson et al. 2013)107 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31)94 LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013) 95 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1) 108 96 ! 109 LOGICAL :: ln_taudif ! logical flag to use the "mean of stress module - module of mean stress" data 110 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 111 REAL(wp) :: rn_efac ! multiplication factor for evaporation 112 REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress 113 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 114 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 115 !!gm ref namelist initialize it so remove the setting to false below 116 LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 117 LOGICAL :: ln_Cd_L15 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 97 LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 98 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 118 99 ! 119 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_atm ! transfer coefficient for momentum (tau) 120 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ch_atm ! transfer coefficient for sensible heat (Q_sens) 121 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ce_atm ! tansfert coefficient for evaporation (Q_lat) 122 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_zu ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 123 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_zu ! air spec. hum. at wind speed height (needed by Lupkes 2015 bulk scheme) 124 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 100 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 101 REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation 102 REAL(wp), PUBLIC :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress 103 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 104 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 105 ! 106 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice 107 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme) 108 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 109 110 LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 111 LOGICAL :: ln_skin_wl ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 112 LOGICAL :: ln_humi_sph ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB 113 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 114 LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB 115 ! 116 INTEGER :: nhumi ! choice of the bulk algorithm 117 ! ! associated indices: 118 INTEGER, PARAMETER :: np_humi_sph = 1 119 INTEGER, PARAMETER :: np_humi_dpt = 2 120 INTEGER, PARAMETER :: np_humi_rlh = 3 125 121 126 122 INTEGER :: nblk ! choice of the bulk algorithm … … 128 124 INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) 129 125 INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003) 130 INTEGER, PARAMETER :: np_COARE_3p 5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013)131 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31)126 INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013) 127 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 45r1) 132 128 133 129 !! * Substitutions … … 144 140 !! *** ROUTINE sbc_blk_alloc *** 145 141 !!------------------------------------------------------------------- 146 ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 147 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 142 ALLOCATE( t_zu(jpi,jpj) , q_zu(jpi,jpj) , & 143 & Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj), & 144 & Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc ) 148 145 ! 149 146 CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) … … 158 155 !! ** Purpose : choose and initialize a bulk formulae formulation 159 156 !! 160 !! ** Method : 157 !! ** Method : 161 158 !! 162 159 !!---------------------------------------------------------------------- 163 INTEGER :: ifpr, jfld! dummy loop indice and argument160 INTEGER :: jfpr ! dummy loop indice and argument 164 161 INTEGER :: ios, ierror, ioptio ! Local integer 165 162 !! 166 163 CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files 167 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i! array of namelist informations on the fields to read164 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 168 165 TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read 169 166 TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " " 170 TYPE(FLD_N) :: sn_slp , sn_ tdif! " "167 TYPE(FLD_N) :: sn_slp , sn_hpgi, sn_hpgj ! " " 171 168 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 172 & sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, & 173 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 174 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 175 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 169 & sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj, & 170 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 171 & cn_dir , rn_zqt, rn_zu, & 172 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, & 173 & ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB 176 174 !!--------------------------------------------------------------------- 177 175 ! … … 179 177 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 180 178 ! 181 ! !** read bulk namelist 179 ! !** read bulk namelist 182 180 REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters 183 181 READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) … … 192 190 ! !** initialization of the chosen bulk formulae (+ check) 193 191 ! !* select the bulk chosen in the namelist and check the choice 194 ioptio = 0 195 IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF 196 IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF 197 IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF 198 IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF 199 ! 192 ioptio = 0 193 IF( ln_NCAR ) THEN 194 nblk = np_NCAR ; ioptio = ioptio + 1 195 ENDIF 196 IF( ln_COARE_3p0 ) THEN 197 nblk = np_COARE_3p0 ; ioptio = ioptio + 1 198 ENDIF 199 IF( ln_COARE_3p6 ) THEN 200 nblk = np_COARE_3p6 ; ioptio = ioptio + 1 201 ENDIF 202 IF( ln_ECMWF ) THEN 203 nblk = np_ECMWF ; ioptio = ioptio + 1 204 ENDIF 200 205 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 206 207 ! !** initialization of the cool-skin / warm-layer parametrization 208 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 209 !! Some namelist sanity tests: 210 IF( ln_NCAR ) & 211 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' ) 212 IF( nn_fsbc /= 1 ) & 213 & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.') 214 END IF 215 216 IF( ln_skin_wl ) THEN 217 !! Check if the frequency of downwelling solar flux input makes sense and if ln_dm2dc=T if it is daily! 218 IF( (sn_qsr%freqh < 0.).OR.(sn_qsr%freqh > 24.) ) & 219 & CALL ctl_stop( 'sbc_blk_init: Warm-layer param. (ln_skin_wl) not compatible with freq. of solar flux > daily' ) 220 IF( (sn_qsr%freqh == 24.).AND.(.NOT. ln_dm2dc) ) & 221 & CALL ctl_stop( 'sbc_blk_init: Please set ln_dm2dc=T for warm-layer param. (ln_skin_wl) to work properly' ) 222 END IF 223 224 ioptio = 0 225 IF( ln_humi_sph ) THEN 226 nhumi = np_humi_sph ; ioptio = ioptio + 1 227 ENDIF 228 IF( ln_humi_dpt ) THEN 229 nhumi = np_humi_dpt ; ioptio = ioptio + 1 230 ENDIF 231 IF( ln_humi_rlh ) THEN 232 nhumi = np_humi_rlh ; ioptio = ioptio + 1 233 ENDIF 234 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 201 235 ! 202 236 IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr 203 237 IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 204 IF( sn_qsr%ln_tint ) THEN 238 IF( sn_qsr%ln_tint ) THEN 205 239 CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', & 206 240 & ' ==> We force time interpolation = .false. for qsr' ) … … 210 244 ! !* set the bulk structure 211 245 ! !- store namelist information in an array 246 IF( ln_blk ) jpfld = 9 247 IF( ln_abl ) jpfld = 11 248 ALLOCATE( slf_i(jpfld) ) 249 ! 212 250 slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj 213 251 slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw 214 252 slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi 215 253 slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow 216 slf_i(jp_slp ) = sn_slp ; slf_i(jp_tdif) = sn_tdif217 !218 lhftau = ln_taudif !- add an extra field if HF stress is used219 jfld = jpfld - COUNT( (/.NOT.lhftau/) )254 slf_i(jp_slp ) = sn_slp 255 IF( ln_abl ) THEN 256 slf_i(jp_hpgi) = sn_hpgi ; slf_i(jp_hpgj) = sn_hpgj 257 END IF 220 258 ! 221 259 ! !- allocate the bulk structure 222 ALLOCATE( sf(j fld), STAT=ierror )260 ALLOCATE( sf(jpfld), STAT=ierror ) 223 261 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 224 DO ifpr= 1, jfld 225 ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 226 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 227 IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 ) & 228 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & 229 & ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 230 262 ! 263 DO jfpr= 1, jpfld 264 ! 265 IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to zero) 266 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 267 sf(jfpr)%fnow(:,:,1) = 0._wp 268 ELSE !-- used field --! 269 IF( ln_abl .AND. & 270 & ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR. & 271 & jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair ) ) THEN ! ABL: some fields are 3D input 272 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) ) 273 IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) ) 274 ELSE ! others or Bulk fields are 2D fiels 275 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 276 IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) 277 ENDIF 278 ! 279 IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 ) & 280 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & 281 & ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 282 ENDIF 231 283 END DO 232 284 ! !- fill the bulk structure with namelist informations 233 285 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 234 286 ! 235 IF 236 !Activated wave module but neither drag nor stokes drift activated237 IF 287 IF( ln_wave ) THEN 288 !Activated wave module but neither drag nor stokes drift activated 289 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN 238 290 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 239 !drag coefficient read from wave model definable only with mfs bulk formulae and core240 ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR ) THEN241 242 ELSEIF 243 291 !drag coefficient read from wave model definable only with mfs bulk formulae and core 292 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN 293 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 294 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN 295 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 244 296 ENDIF 245 297 ELSE 246 IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 247 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 248 & 'with drag coefficient (ln_cdgw =T) ' , & 249 & 'or Stokes Drift (ln_sdw=T) ' , & 250 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', & 251 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 252 ENDIF 253 ! 254 ! 298 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 299 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 300 & 'with drag coefficient (ln_cdgw =T) ' , & 301 & 'or Stokes Drift (ln_sdw=T) ' , & 302 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', & 303 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 304 ENDIF 305 ! 306 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient 307 rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level 308 rn_zu = ght_abl(2) 309 IF(lwp) WRITE(numout,*) 310 IF(lwp) WRITE(numout,*) ' ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude' 311 ENDIF 312 ! 313 ! set transfer coefficients to default sea-ice values 314 Cd_ice(:,:) = rCd_ice 315 Ch_ice(:,:) = rCd_ice 316 Ce_ice(:,:) = rCd_ice 317 ! 255 318 IF(lwp) THEN !** Control print 256 319 ! 257 WRITE(numout,*) !* namelist 320 WRITE(numout,*) !* namelist 258 321 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 259 322 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 260 323 WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0 261 WRITE(numout,*) ' "COARE 3.5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p0 262 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31) ln_ECMWF = ', ln_ECMWF 263 WRITE(numout,*) ' add High freq.contribution to the stress module ln_taudif = ', ln_taudif 324 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 325 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 264 326 WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt 265 327 WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu … … 275 337 CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)' 276 338 CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)' 277 CASE( np_COARE_3p 5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)'278 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)'339 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 340 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)' 279 341 END SELECT 280 342 ! 343 WRITE(numout,*) 344 WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs 345 WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl 346 ! 347 WRITE(numout,*) 348 SELECT CASE( nhumi ) !* Print the choice of air humidity 349 CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]' 350 CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]' 351 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 352 END SELECT 353 ! 281 354 ENDIF 282 355 ! … … 291 364 !! (momentum, heat, freshwater and runoff) 292 365 !! 293 !! ** Method : (1) READ each fluxes in NetCDF files: 294 !! the 10m wind velocity (i-component) (m/s) at T-point 295 !! the 10m wind velocity (j-component) (m/s) at T-point 296 !! the 10m or 2m specific humidity ( % ) 297 !! the solar heat (W/m2) 298 !! the Long wave (W/m2) 299 !! the 10m or 2m air temperature (Kelvin) 300 !! the total precipitation (rain+snow) (Kg/m2/s) 301 !! the snow (solid prcipitation) (kg/m2/s) 302 !! the tau diff associated to HF tau (N/m2) at T-point (ln_taudif=T) 303 !! (2) CALL blk_oce 366 !! ** Method : 367 !! (1) READ each fluxes in NetCDF files: 368 !! the wind velocity (i-component) at z=rn_zu (m/s) at T-point 369 !! the wind velocity (j-component) at z=rn_zu (m/s) at T-point 370 !! the specific humidity at z=rn_zqt (kg/kg) 371 !! the air temperature at z=rn_zqt (Kelvin) 372 !! the solar heat (W/m2) 373 !! the Long wave (W/m2) 374 !! the total precipitation (rain+snow) (Kg/m2/s) 375 !! the snow (solid precipitation) (kg/m2/s) 376 !! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) 377 !! (2) CALL blk_oce_1 and blk_oce_2 304 378 !! 305 379 !! C A U T I O N : never mask the surface stress fields … … 318 392 !!---------------------------------------------------------------------- 319 393 INTEGER, INTENT(in) :: kt ! ocean time step 320 !!--------------------------------------------------------------------- 394 !!---------------------------------------------------------------------- 395 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zevp 396 REAL(wp) :: ztmp 397 !!---------------------------------------------------------------------- 321 398 ! 322 399 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 323 ! 400 401 ! Sanity/consistence test on humidity at first time step to detect potential screw-up: 402 IF( kt == nit000 ) THEN 403 WRITE(numout,*) '' 404 #if defined key_agrif 405 WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 406 #else 407 ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain 408 IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points! 409 ztmp = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztmp ! mean humidity over ocean on proc 410 SELECT CASE( nhumi ) 411 CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!) 412 IF( (ztmp < 0._wp) .OR. (ztmp > 0.065) ) ztmp = -1._wp 413 CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K] 414 IF( (ztmp < 110._wp).OR.(ztmp > 320._wp) ) ztmp = -1._wp 415 CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%] 416 IF( (ztmp < 0._wp) .OR.(ztmp > 100._wp) ) ztmp = -1._wp 417 END SELECT 418 IF(ztmp < 0._wp) THEN 419 WRITE(numout,'(" Mean humidity value found on proc #",i5.5," is: ",f)'), narea, ztmp 420 CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', & 421 & ' ==> check the unit in your input files' , & 422 & ' ==> check consistence of namelist choice: specific? relative? dew-point?', & 423 & ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' ) 424 END IF 425 END IF 426 WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 427 #endif 428 WRITE(numout,*) '' 429 END IF !IF( kt == nit000 ) 324 430 ! ! compute the surface ocean fluxes using bulk formulea 325 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 326 431 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 432 CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! <<= in 433 & sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & ! <<= in 434 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 435 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 436 & zssq, zcd_du, zsen, zevp ) ! =>> out 437 438 CALL blk_oce_2( sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in 439 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), & ! <<= in 440 & sf(jp_snow)%fnow(:,:,1), sst_m, & ! <<= in 441 & zsen, zevp ) ! <=> in out 442 ENDIF 443 ! 327 444 #if defined key_cice 328 445 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 329 446 qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) 330 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 331 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 332 ENDIF 447 IF( ln_dm2dc ) THEN 448 qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 449 ELSE 450 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 451 ENDIF 333 452 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 334 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 453 454 SELECT CASE( nhumi ) 455 CASE( np_humi_sph ) 456 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 457 CASE( np_humi_dpt ) 458 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 459 CASE( np_humi_rlh ) 460 qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file 461 END SELECT 462 335 463 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 336 464 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac … … 343 471 344 472 345 SUBROUTINE blk_oce( kt, sf, pst, pu, pv ) 346 !!--------------------------------------------------------------------- 347 !! *** ROUTINE blk_oce *** 348 !! 349 !! ** Purpose : provide the momentum, heat and freshwater fluxes at 350 !! the ocean surface at each time step 351 !! 352 !! ** Method : bulk formulea for the ocean using atmospheric 353 !! fields read in sbc_read 473 SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, & ! inp 474 & pslp , pst , pu , pv, & ! inp 475 & pqsr , pqlw , & ! inp 476 & pssq , pcd_du, psen , pevp ) ! out 477 !!--------------------------------------------------------------------- 478 !! *** ROUTINE blk_oce_1 *** 479 !! 480 !! ** Purpose : if ln_blk=T, computes surface momentum, heat and freshwater fluxes 481 !! if ln_abl=T, computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 482 !! 483 !! ** Method : bulk formulae using atmospheric fields from : 484 !! if ln_blk=T, atmospheric fields read in sbc_read 485 !! if ln_abl=T, the ABL model at previous time-step 486 !! 487 !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) 488 !! - pcd_du : Cd x |dU| at T-points (m/s) 489 !! - psen : Ch x |dU| at T-points (m/s) 490 !! - pevp : Ce x |dU| at T-points (m/s) 491 !!--------------------------------------------------------------------- 492 INTEGER , INTENT(in ) :: kt ! time step index 493 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s] 494 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s] 495 REAL(wp), INTENT(in ), DIMENSION(:,:) :: phumi ! specific humidity at T-points [kg/kg] 496 REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin] 497 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa] 498 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celcius] 499 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 500 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 501 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr ! 502 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw ! 503 REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] 504 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] 505 REAL(wp), INTENT( out), DIMENSION(:,:) :: psen ! Ch x |dU| at T-points [m/s] 506 REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp ! Ce x |dU| at T-points [m/s] 507 ! 508 INTEGER :: ji, jj ! dummy loop indices 509 REAL(wp) :: zztmp ! local variable 510 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 511 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 512 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 513 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] 514 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] 515 REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean 516 REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean 517 REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean 518 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux 519 REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2 520 !!--------------------------------------------------------------------- 521 ! 522 ! local scalars ( place there for vector optimisation purposes) 523 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 524 525 ! ----------------------------------------------------------------------------- ! 526 ! 0 Wind components and module at T-point relative to the moving ocean ! 527 ! ----------------------------------------------------------------------------- ! 528 529 ! ... components ( U10m - U_oce ) at T-point (unmasked) 530 #if defined key_cyclone 531 zwnd_i(:,:) = 0._wp 532 zwnd_j(:,:) = 0._wp 533 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 534 DO jj = 2, jpjm1 535 DO ji = fs_2, fs_jpim1 ! vect. opt. 536 pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 537 pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 538 END DO 539 END DO 540 #endif 541 DO jj = 2, jpjm1 542 DO ji = fs_2, fs_jpim1 ! vect. opt. 543 zwnd_i(ji,jj) = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 544 zwnd_j(ji,jj) = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 545 END DO 546 END DO 547 CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 548 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 549 wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) & 550 & + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1) 551 552 ! ----------------------------------------------------------------------------- ! 553 ! I Solar FLUX ! 554 ! ----------------------------------------------------------------------------- ! 555 556 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 557 zztmp = 1. - albo 558 IF( ln_dm2dc ) THEN 559 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 560 ELSE 561 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 562 ENDIF 563 564 565 ! ----------------------------------------------------------------------------- ! 566 ! II Turbulent FLUXES ! 567 ! ----------------------------------------------------------------------------- ! 568 569 ! specific humidity at SST 570 pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) ) 571 572 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 573 zztmp1(:,:) = zst(:,:) 574 zztmp2(:,:) = pssq(:,:) 575 ENDIF 576 577 ! specific humidity of air at "rn_zqt" m above the sea 578 SELECT CASE( nhumi ) 579 CASE( np_humi_sph ) 580 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 581 CASE( np_humi_dpt ) 582 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 583 zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) ) 584 CASE( np_humi_rlh ) 585 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 586 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 587 END SELECT 588 ! 589 ! potential temperature of air at "rn_zqt" m above the sea 590 IF( ln_abl ) THEN 591 ztpot = ptair(:,:) 592 ELSE 593 ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 594 ! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 595 ! (since reanalysis products provide T at z, not theta !) 596 !#LB: because AGRIF hates functions that return something else than a scalar, need to 597 ! use scalar version of gamma_moist() ... 598 DO jj = 1, jpj 599 DO ji = 1, jpi 600 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 601 END DO 602 END DO 603 ENDIF 604 605 606 607 !! Time to call the user-selected bulk parameterization for 608 !! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more... 609 SELECT CASE( nblk ) 610 611 CASE( np_NCAR ) 612 CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, & 613 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 614 615 CASE( np_COARE_3p0 ) 616 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 617 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 618 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 619 620 CASE( np_COARE_3p6 ) 621 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 622 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 623 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 624 625 CASE( np_ECMWF ) 626 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 627 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 628 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 629 630 CASE DEFAULT 631 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 632 633 END SELECT 634 635 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 636 !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq: 637 WHERE ( fr_i < 0.001_wp ) 638 ! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!! 639 zst(:,:) = zst(:,:)*tmask(:,:,1) 640 pssq(:,:) = pssq(:,:)*tmask(:,:,1) 641 ELSEWHERE 642 ! we forget about the update... 643 zst(:,:) = zztmp1(:,:) !#LB: using what we backed up before skin-algo 644 pssq(:,:) = zztmp2(:,:) !#LB: " " " 645 END WHERE 646 END IF 647 648 !! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef. 649 !! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef. 650 651 IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 652 !! If zu == zt, then ensuring once for all that: 653 t_zu(:,:) = ztpot(:,:) 654 q_zu(:,:) = zqair(:,:) 655 ENDIF 656 657 658 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90 659 ! ------------------------------------------------------------- 660 661 IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 662 !! FL do we need this multiplication by tmask ... ??? 663 DO jj = 1, jpj 664 DO ji = 1, jpi 665 zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 666 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 667 pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 668 psen(ji,jj) = zztmp * zch_oce(ji,jj) 669 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 670 END DO 671 END DO 672 ELSE !== BLK formulation ==! turbulent fluxes computation 673 CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 674 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 675 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 676 & taum(:,:), psen(:,:), zqla(:,:), & 677 & pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 678 679 zqla(:,:) = zqla(:,:) * tmask(:,:,1) 680 psen(:,:) = psen(:,:) * tmask(:,:,1) 681 taum(:,:) = taum(:,:) * tmask(:,:,1) 682 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 683 684 ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 685 zcd_oce = 0._wp 686 WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 687 zwnd_i = zcd_oce * zwnd_i 688 zwnd_j = zcd_oce * zwnd_j 689 690 CALL iom_put( "taum_oce", taum ) ! output wind stress module 691 692 ! ... utau, vtau at U- and V_points, resp. 693 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 694 ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 695 DO jj = 1, jpjm1 696 DO ji = 1, fs_jpim1 697 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & 698 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 699 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) & 700 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 701 END DO 702 END DO 703 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 704 705 IF(ln_ctl) THEN 706 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') 707 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & 708 & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 709 ENDIF 710 ! 711 ENDIF 712 ! 713 IF(ln_ctl) THEN 714 CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' ) 715 CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' ) 716 CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' ) 717 ENDIF 718 ! 719 END SUBROUTINE blk_oce_1 720 721 722 SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec, & ! <<= in 723 & psnow, pst , psen, pevp ) ! <<= in 724 !!--------------------------------------------------------------------- 725 !! *** ROUTINE blk_oce_2 *** 726 !! 727 !! ** Purpose : finalize the momentum, heat and freshwater fluxes computation 728 !! at the ocean surface at each time step knowing Cd, Ch, Ce and 729 !! atmospheric variables (from ABL or external data) 354 730 !! 355 731 !! ** Outputs : - utau : i-component of the stress at U-point (N/m2) … … 360 736 !! - qns : Non Solar heat flux over the ocean (W/m2) 361 737 !! - emp : evaporation minus precipitation (kg/m2/s) 362 !! 363 !! ** Nota : sf has to be a dummy argument for AGRIF on NEC 364 !!--------------------------------------------------------------------- 365 INTEGER , INTENT(in ) :: kt ! time step index 366 TYPE(fld), INTENT(inout), DIMENSION(:) :: sf ! input data 367 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pst ! surface temperature [Celcius] 368 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 369 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 738 !!--------------------------------------------------------------------- 739 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair 740 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr 741 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw 742 REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec 743 REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow 744 REAL(wp), INTENT(in), DIMENSION(:,:) :: pst ! surface temperature [Celcius] 745 REAL(wp), INTENT(in), DIMENSION(:,:) :: psen 746 REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp 370 747 ! 371 748 INTEGER :: ji, jj ! dummy loop indices 372 REAL(wp) :: zztmp ! local variable 373 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 374 REAL(wp), DIMENSION(jpi,jpj) :: zsq ! specific humidity at pst 375 REAL(wp), DIMENSION(jpi,jpj) :: zqlw, zqsb ! long wave and sensible heat fluxes 376 REAL(wp), DIMENSION(jpi,jpj) :: zqla, zevap ! latent heat fluxes and evaporation 749 REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable 750 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes 751 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation 377 752 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 378 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]379 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K]380 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! density of air [kg/m^3]381 753 !!--------------------------------------------------------------------- 382 754 ! … … 384 756 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 385 757 758 386 759 ! ----------------------------------------------------------------------------- ! 387 ! 0 Wind components and module at T-point relative to the moving ocean!760 ! III Net longwave radiative FLUX ! 388 761 ! ----------------------------------------------------------------------------- ! 389 762 390 ! ... components ( U10m - U_oce ) at T-point (unmasked) 391 !!gm move zwnd_i (_j) set to zero inside the key_cyclone ??? 392 zwnd_i(:,:) = 0._wp 393 zwnd_j(:,:) = 0._wp 394 #if defined key_cyclone 395 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 396 DO jj = 2, jpjm1 397 DO ji = fs_2, fs_jpim1 ! vect. opt. 398 sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj) 399 sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj) 400 END DO 401 END DO 402 #endif 403 DO jj = 2, jpjm1 404 DO ji = fs_2, fs_jpim1 ! vect. opt. 405 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 406 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 407 END DO 408 END DO 409 CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 410 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 411 wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) & 412 & + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1) 413 414 ! ----------------------------------------------------------------------------- ! 415 ! I Radiative FLUXES ! 416 ! ----------------------------------------------------------------------------- ! 417 418 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 419 zztmp = 1. - albo 420 IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 421 ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 422 ENDIF 423 424 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 425 426 ! ----------------------------------------------------------------------------- ! 427 ! II Turbulent FLUXES ! 428 ! ----------------------------------------------------------------------------- ! 429 430 ! ... specific humidity at SST and IST tmask( 431 zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 432 !! 433 !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 434 !! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 435 !! (since reanalysis products provide T at z, not theta !) 436 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 437 438 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 439 ! 440 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2 441 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 442 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0 443 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 444 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5 445 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 446 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 447 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 448 CASE DEFAULT 449 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 450 END SELECT 451 452 ! ! Compute true air density : 453 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 454 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 455 ELSE ! At zt: 456 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 457 END IF 458 459 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. 460 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 461 462 DO jj = 1, jpj ! tau module, i and j component 463 DO ji = 1, jpi 464 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 465 taum (ji,jj) = zztmp * wndm (ji,jj) 466 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 467 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 468 END DO 469 END DO 470 471 ! ! add the HF tau contribution to the wind stress module 472 IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 473 474 CALL iom_put( "taum_oce", taum ) ! output wind stress module 475 476 ! ... utau, vtau at U- and V_points, resp. 477 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 478 ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 479 DO jj = 1, jpjm1 480 DO ji = 1, fs_jpim1 481 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & 482 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 483 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) & 484 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 485 END DO 486 END DO 487 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 763 !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 764 !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 765 zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux 488 766 489 767 ! Turbulent fluxes over ocean 490 768 ! ----------------------------- 491 769 492 ! zqla used as temporary array, for rho*U (common term of bulk formulae): 493 zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1) 494 495 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 496 !! q_air and t_air are given at 10m (wind reference height) 497 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 498 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 499 ELSE 500 !! q_air and t_air are not given at 10m (wind reference height) 501 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 502 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 503 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed 504 ENDIF 505 506 zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux 507 770 ! use scalar version of L_vap() for AGRIF compatibility 771 DO jj = 1, jpj 772 DO ji = 1, jpi 773 zqla(ji,jj) = L_vap( zst(ji,jj) ) * pevp(ji,jj) * -1._wp ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 774 ENDDO 775 ENDDO 508 776 509 777 IF(ln_ctl) THEN 510 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' ) 511 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' ) 512 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 513 CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' ) 514 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, & 515 & tab2d_2=vtau , clinfo2= ' vtau : ', mask2=vmask ) 516 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce: wndm : ') 517 CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce: zst : ') 778 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' ) 779 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 780 518 781 ENDIF 519 782 520 783 ! ----------------------------------------------------------------------------- ! 521 ! I IITotal FLUXES !784 ! IV Total FLUXES ! 522 785 ! ----------------------------------------------------------------------------- ! 523 786 ! 524 emp (:,:) = ( zevap(:,:)& ! mass flux (evap. - precip.)525 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1)526 ! 527 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)& ! Downward Non Solar528 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus & ! remove latent melting heat for solid precip529 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST530 & + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac& ! add liquid precip heat content at Tair531 & * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp &532 & + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)533 & * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi787 emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.) 788 & - pprec(:,:) * rn_pfac ) * tmask(:,:,1) 789 ! 790 qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar 791 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 792 & - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST !LB??? pst is Celsius !? 793 & + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair 794 & * ( ptair(:,:) - rt0 ) * rcp & 795 & + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 796 & * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 534 797 qns(:,:) = qns(:,:) * tmask(:,:,1) 535 798 ! 536 799 #if defined key_si3 537 qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)! non solar without emp (only needed by SI3)800 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3) 538 801 qsr_oce(:,:) = qsr(:,:) 539 802 #endif 540 803 ! 804 CALL iom_put( "rho_air" , rhoa*tmask(:,:,1) ) ! output air density [kg/m^3] 805 CALL iom_put( "evap_oce" , pevp ) ! evaporation 806 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 807 CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean 808 CALL iom_put( "qla_oce" , zqla ) ! output downward latent heat over the ocean 809 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 810 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 811 CALL iom_put( 'snowpre', sprecip ) ! Snow 812 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 813 ! 541 814 IF ( nn_ice == 0 ) THEN 542 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 543 CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean 544 CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean 545 CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean 546 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 547 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 548 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 549 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 550 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 551 CALL iom_put( 'snowpre', sprecip ) ! Snow 552 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 815 CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla ) ! output downward heat content of E-P over the ocean 816 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 817 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 818 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 819 ENDIF 820 ! 821 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 822 CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius 823 CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference... 553 824 ENDIF 554 825 ! 555 826 IF(ln_ctl) THEN 556 CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=zqlw , clinfo2=' zqlw : ') 557 CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') 558 CALL prt_ctl(tab2d_1=pst , clinfo1=' blk_oce: pst : ', tab2d_2=emp , clinfo2=' emp : ') 559 CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, & 560 & tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask ) 561 ENDIF 562 ! 563 END SUBROUTINE blk_oce 564 565 566 567 FUNCTION rho_air( ptak, pqa, pslp ) 568 !!------------------------------------------------------------------------------- 569 !! *** FUNCTION rho_air *** 570 !! 571 !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 572 !! 573 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 574 !!------------------------------------------------------------------------------- 575 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 576 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 577 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa] 578 REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3] 579 !!------------------------------------------------------------------------------- 580 ! 581 rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) ) 582 ! 583 END FUNCTION rho_air 584 585 586 FUNCTION cp_air( pqa ) 587 !!------------------------------------------------------------------------------- 588 !! *** FUNCTION cp_air *** 589 !! 590 !! ** Purpose : Compute specific heat (Cp) of moist air 591 !! 592 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 593 !!------------------------------------------------------------------------------- 594 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 595 REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg] 596 !!------------------------------------------------------------------------------- 597 ! 598 Cp_air = Cp_dry + Cp_vap * pqa 599 ! 600 END FUNCTION cp_air 601 602 603 FUNCTION q_sat( ptak, pslp ) 604 !!---------------------------------------------------------------------------------- 605 !! *** FUNCTION q_sat *** 606 !! 607 !! ** Purpose : Specific humidity at saturation in [kg/kg] 608 !! Based on accurate estimate of "e_sat" 609 !! aka saturation water vapor (Goff, 1957) 610 !! 611 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 612 !!---------------------------------------------------------------------------------- 613 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 614 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa] 615 REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg] 616 ! 617 INTEGER :: ji, jj ! dummy loop indices 618 REAL(wp) :: ze_sat, ztmp ! local scalar 619 !!---------------------------------------------------------------------------------- 620 ! 621 DO jj = 1, jpj 622 DO ji = 1, jpi 623 ! 624 ztmp = rt0 / ptak(ji,jj) 625 ! 626 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 627 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) & 628 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) & 629 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 ) 630 ! 631 q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat ) ! 0.01 because SLP is in [Pa] 632 ! 633 END DO 634 END DO 635 ! 636 END FUNCTION q_sat 637 638 639 FUNCTION gamma_moist( ptak, pqa ) 640 !!---------------------------------------------------------------------------------- 641 !! *** FUNCTION gamma_moist *** 642 !! 643 !! ** Purpose : Compute the moist adiabatic lapse-rate. 644 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 645 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 646 !! 647 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 648 !!---------------------------------------------------------------------------------- 649 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 650 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg] 651 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate 652 ! 653 INTEGER :: ji, jj ! dummy loop indices 654 REAL(wp) :: zrv, ziRT ! local scalar 655 !!---------------------------------------------------------------------------------- 656 ! 657 DO jj = 1, jpj 658 DO ji = 1, jpi 659 zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 660 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT 661 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 662 END DO 663 END DO 664 ! 665 END FUNCTION gamma_moist 666 667 668 FUNCTION L_vap( psst ) 669 !!--------------------------------------------------------------------------------- 670 !! *** FUNCTION L_vap *** 671 !! 672 !! ** Purpose : Compute the latent heat of vaporization of water from temperature 673 !! 674 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 675 !!---------------------------------------------------------------------------------- 676 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg] 677 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K] 678 !!---------------------------------------------------------------------------------- 679 ! 680 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6 681 ! 682 END FUNCTION L_vap 827 CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ') 828 CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') 829 CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ') 830 ENDIF 831 ! 832 END SUBROUTINE blk_oce_2 833 683 834 684 835 #if defined key_si3 … … 686 837 !! 'key_si3' SI3 sea-ice model 687 838 !!---------------------------------------------------------------------- 688 !! blk_ice_ tau: provide the air-ice stress689 !! blk_ice_ flx: provide the heat and mass fluxes at air-ice interface839 !! blk_ice_1 : provide the air-ice stress 840 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 690 841 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 691 842 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 692 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 843 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 693 844 !!---------------------------------------------------------------------- 694 845 695 SUBROUTINE blk_ice_tau 696 !!--------------------------------------------------------------------- 697 !! *** ROUTINE blk_ice_tau *** 846 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui, & ! inputs 847 & putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs 848 !!--------------------------------------------------------------------- 849 !! *** ROUTINE blk_ice_1 *** 698 850 !! 699 851 !! ** Purpose : provide the surface boundary condition over sea-ice … … 703 855 !! NB: ice drag coefficient is assumed to be a constant 704 856 !!--------------------------------------------------------------------- 857 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa] 858 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s] 859 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s] 860 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s] 861 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: phumi ! atmospheric wind at T-point [m/s] 862 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s] 863 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! " 864 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K] 865 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk 866 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk 867 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl 868 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl 869 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl 870 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl 871 ! 705 872 INTEGER :: ji, jj ! dummy loop indices 706 REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point707 873 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 708 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau) 709 !!--------------------------------------------------------------------- 710 ! 711 ! set transfer coefficients to default sea-ice values 712 Cd_atm(:,:) = Cd_ice 713 Ch_atm(:,:) = Cd_ice 714 Ce_atm(:,:) = Cd_ice 715 716 wndm_ice(:,:) = 0._wp !!gm brutal.... 874 REAL(wp) :: zootm_su ! sea-ice surface mean temperature 875 REAL(wp) :: zztmp1, zztmp2 ! temporary arrays 876 REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau) 877 !!--------------------------------------------------------------------- 878 ! 717 879 718 880 ! ------------------------------------------------------------ ! … … 722 884 DO jj = 2, jpjm1 723 885 DO ji = fs_2, fs_jpim1 ! vect. opt. 724 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )725 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )886 zwndi_t = ( pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj ) + puice(ji,jj) ) ) 887 zwndj_t = ( pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji ,jj-1) + pvice(ji,jj) ) ) 726 888 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 727 889 END DO … … 731 893 ! Make ice-atm. drag dependent on ice concentration 732 894 IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations 733 CALL Cdn10_Lupkes2012( Cd_atm ) 734 Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical 895 CALL Cdn10_Lupkes2012( Cd_ice ) 896 Ch_ice(:,:) = Cd_ice(:,:) ! momentum and heat transfer coef. are considered identical 897 Ce_ice(:,:) = Cd_ice(:,:) 735 898 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 736 CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 737 ENDIF 738 739 !! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef. 740 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef. 899 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 900 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 901 ENDIF 902 903 !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef. 904 !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef. 741 905 742 906 ! local scalars ( place there for vector optimisation purposes) 743 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 744 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 745 746 !!gm brutal.... 747 utau_ice (:,:) = 0._wp 748 vtau_ice (:,:) = 0._wp 749 !!gm end 750 751 ! ------------------------------------------------------------ ! 752 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 753 ! ------------------------------------------------------------ ! 754 ! C-grid ice dynamics : U & V-points (same as ocean) 755 DO jj = 2, jpjm1 756 DO ji = fs_2, fs_jpim1 ! vect. opt. 757 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 758 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 759 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 760 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 907 !IF (ln_abl) rhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI) 908 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 909 910 IF( ln_blk ) THEN 911 ! ------------------------------------------------------------ ! 912 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 913 ! ------------------------------------------------------------ ! 914 ! C-grid ice dynamics : U & V-points (same as ocean) 915 DO jj = 2, jpjm1 916 DO ji = fs_2, fs_jpim1 ! vect. opt. 917 putaui(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) & 918 & + rhoa(ji ,jj) * zcd_dui(ji ,jj) ) & 919 & * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 920 pvtaui(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) & 921 & + rhoa(ji,jj ) * zcd_dui(ji,jj ) ) & 922 & * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 923 END DO 761 924 END DO 762 END DO 763 CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. ) 764 ! 765 ! 766 IF(ln_ctl) THEN 767 CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ') 768 CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') 769 ENDIF 770 ! 771 END SUBROUTINE blk_ice_tau 772 773 774 SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 775 !!--------------------------------------------------------------------- 776 !! *** ROUTINE blk_ice_flx *** 925 CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 926 ! 927 IF(ln_ctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' & 928 & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) 929 ELSE 930 zztmp1 = 11637800.0_wp 931 zztmp2 = -5897.8_wp 932 DO jj = 1, jpj 933 DO ji = 1, jpi 934 pcd_dui(ji,jj) = zcd_dui (ji,jj) 935 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 936 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 937 zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??) 938 pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 939 END DO 940 END DO 941 ENDIF 942 ! 943 IF(ln_ctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') 944 ! 945 END SUBROUTINE blk_ice_1 946 947 948 SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow ) 949 !!--------------------------------------------------------------------- 950 !! *** ROUTINE blk_ice_2 *** 777 951 !! 778 952 !! ** Purpose : provide the heat and mass fluxes at air-ice interface … … 784 958 !! caution : the net upward water flux has with mm/day unit 785 959 !!--------------------------------------------------------------------- 786 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature 960 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K] 787 961 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness 788 962 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 789 963 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 964 REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair 965 REAL(wp), DIMENSION(:,: ), INTENT(in) :: phumi 966 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp 967 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqlw 968 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec 969 REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow 790 970 !! 791 971 INTEGER :: ji, jj, jl ! dummy loop indices 792 972 REAL(wp) :: zst3 ! local variable 793 973 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 794 REAL(wp) :: zztmp, z 1_rLsub! - -974 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 795 975 REAL(wp) :: zfr1, zfr2 ! local variables 796 976 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature … … 800 980 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 801 981 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3) 802 REAL(wp), DIMENSION(jpi,jpj) :: z rhoa982 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB 803 983 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2 804 984 !!--------------------------------------------------------------------- 805 985 ! 806 zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars 807 zcoef_dqla = -Ls * 11637800. * (-5897.8) 808 ! 809 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 986 zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars 987 zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 988 ! 989 SELECT CASE( nhumi ) 990 CASE( np_humi_sph ) 991 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 992 CASE( np_humi_dpt ) 993 zqair(:,:) = q_sat( phumi(:,:), pslp ) 994 CASE( np_humi_rlh ) 995 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 996 END SELECT 810 997 ! 811 998 zztmp = 1. / ( 1. - albo ) 812 WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 813 ELSEWHERE ; z1_st(:,:,:) = 0._wp 999 WHERE( ptsu(:,:,:) /= 0._wp ) 1000 z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 1001 ELSEWHERE 1002 z1_st(:,:,:) = 0._wp 814 1003 END WHERE 815 1004 ! ! ========================== ! … … 825 1014 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 826 1015 ! Long Wave (lw) 827 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef* ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)1016 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 828 1017 ! lw sensitivity 829 1018 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 … … 833 1022 ! ----------------------------! 834 1023 835 ! ... turbulent heat fluxes with Ch_ atm recalculated in blk_ice_tau1024 ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 836 1025 ! Sensible Heat 837 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))1026 z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj)) 838 1027 ! Latent Heat 839 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 840 & ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 1028 zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 1029 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1030 & ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 841 1031 ! Latent heat sensitivity for ice (Dqla/Dt) 842 1032 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 843 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * C h_atm(ji,jj) * wndm_ice(ji,jj) * &844 & z1_st(ji,jj,jl) *z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))1033 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1034 & z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2 845 1035 ELSE 846 1036 dqla_ice(ji,jj,jl) = 0._wp … … 848 1038 849 1039 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 850 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)1040 z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 851 1041 852 1042 ! ----------------------------! … … 863 1053 END DO 864 1054 ! 865 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]866 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]867 CALL iom_put( 'snowpre', sprecip ) 868 CALL iom_put( 'precip' , tprecip ) 1055 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s] 1056 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s] 1057 CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation 1058 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 869 1059 870 1060 ! --- evaporation --- ! … … 883 1073 ! --- heat flux associated with emp --- ! 884 1074 qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 885 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp& ! liquid precip at Tair1075 & + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair 886 1076 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) 887 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1077 & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 888 1078 qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only) 889 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1079 & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 890 1080 891 1081 ! --- total solar and non solar fluxes --- ! … … 895 1085 896 1086 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 897 qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1087 qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 898 1088 899 1089 ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 900 1090 DO jl = 1, jpl 901 1091 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 902 ! ! But we do not have Tice => consider it at 0degC => evap=0 1092 ! ! But we do not have Tice => consider it at 0degC => evap=0 903 1093 END DO 904 1094 … … 907 1097 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 908 1098 ! 909 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1099 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 910 1100 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 911 1101 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 912 1102 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 913 1103 ELSEWHERE ! zero when hs>0 914 qtr_ice_top(:,:,:) = 0._wp 1104 qtr_ice_top(:,:,:) = 0._wp 915 1105 END WHERE 916 1106 ! … … 944 1134 ENDIF 945 1135 ! 946 END SUBROUTINE blk_ice_ flx947 1136 END SUBROUTINE blk_ice_2 1137 948 1138 949 1139 SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) … … 954 1144 !! to force sea ice / snow thermodynamics 955 1145 !! in the case conduction flux is emulated 956 !! 1146 !! 957 1147 !! ** Method : compute surface energy balance assuming neglecting heat storage 958 1148 !! following the 0-layer Semtner (1976) approach … … 979 1169 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor 980 1170 !!--------------------------------------------------------------------- 981 1171 982 1172 ! -------------------------------------! 983 1173 ! I Enhanced conduction factor ! … … 987 1177 ! 988 1178 zgfac(:,:,:) = 1._wp 989 1179 990 1180 IF( ld_virtual_itd ) THEN 991 1181 ! … … 993 1183 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 994 1184 zfac3 = 2._wp / zepsilon 995 ! 996 DO jl = 1, jpl 1185 ! 1186 DO jl = 1, jpl 997 1187 DO jj = 1 , jpj 998 1188 DO ji = 1, jpi … … 1002 1192 END DO 1003 1193 END DO 1004 ! 1005 ENDIF 1006 1194 ! 1195 ENDIF 1196 1007 1197 ! -------------------------------------------------------------! 1008 1198 ! II Surface temperature and conduction flux ! … … 1014 1204 DO jj = 1 , jpj 1015 1205 DO ji = 1, jpi 1016 ! 1206 ! 1017 1207 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1018 1208 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) … … 1031 1221 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1032 1222 qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) ) & 1033 1223 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1034 1224 1035 1225 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1036 hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 1226 hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 1037 1227 1038 1228 END DO 1039 1229 END DO 1040 1230 ! 1041 END DO 1042 ! 1231 END DO 1232 ! 1043 1233 END SUBROUTINE blk_ice_qcn 1044 1045 1046 SUBROUTINE Cdn10_Lupkes2012( Cd )1234 1235 1236 SUBROUTINE Cdn10_Lupkes2012( pcd ) 1047 1237 !!---------------------------------------------------------------------- 1048 1238 !! *** ROUTINE Cdn10_Lupkes2012 *** 1049 1239 !! 1050 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1240 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1051 1241 !! to make it dependent on edges at leads, melt ponds and flows. 1052 1242 !! After some approximations, this can be resumed to a dependency 1053 1243 !! on ice concentration. 1054 !! 1244 !! 1055 1245 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50) 1056 1246 !! with the highest level of approximation: level4, eq.(59) … … 1064 1254 !! 1065 1255 !! This new drag has a parabolic shape (as a function of A) starting at 1066 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1256 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1067 1257 !! and going down to Cdi(say 1.4e-3) for A=1 1068 1258 !! … … 1074 1264 !! 1075 1265 !!---------------------------------------------------------------------- 1076 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd1266 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd 1077 1267 REAL(wp), PARAMETER :: zCe = 2.23e-03_wp 1078 1268 REAL(wp), PARAMETER :: znu = 1._wp … … 1089 1279 1090 1280 ! ice-atm drag 1091 Cd(:,:) =Cd_ice + & ! pure ice drag1092 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology1093 1281 pcd(:,:) = rCd_ice + & ! pure ice drag 1282 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1283 1094 1284 END SUBROUTINE Cdn10_Lupkes2012 1095 1285 1096 1286 1097 SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )1287 SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch ) 1098 1288 !!---------------------------------------------------------------------- 1099 1289 !! *** ROUTINE Cdn10_Lupkes2015 *** 1100 1290 !! 1101 1291 !! ** pUrpose : Alternative turbulent transfert coefficients formulation 1102 !! between sea-ice and atmosphere with distinct momentum 1103 !! and heat coefficients depending on sea-ice concentration 1292 !! between sea-ice and atmosphere with distinct momentum 1293 !! and heat coefficients depending on sea-ice concentration 1104 1294 !! and atmospheric stability (no meltponds effect for now). 1105 !! 1295 !! 1106 1296 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1107 1297 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1108 1298 !! it considers specific skin and form drags (Andreas et al. 2010) 1109 !! to compute neutral transfert coefficients for both heat and 1299 !! to compute neutral transfert coefficients for both heat and 1110 1300 !! momemtum fluxes. Atmospheric stability effect on transfert 1111 1301 !! coefficient is also taken into account following Louis (1979). … … 1116 1306 !!---------------------------------------------------------------------- 1117 1307 ! 1118 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd 1119 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch 1120 REAL(wp), DIMENSION(jpi,jpj) :: ztm_su, zst, zqo_sat, zqi_sat 1308 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K] 1309 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa] 1310 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient 1311 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient 1312 REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat 1121 1313 ! 1122 1314 ! ECHAM6 constants … … 1146 1338 !!---------------------------------------------------------------------- 1147 1339 1148 ! mean temperature1149 WHERE( at_i_b(:,:) > 1.e-20 ) ; ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)1150 ELSEWHERE ; ztm_su(:,:) = rt01151 ENDWHERE1152 1153 1340 ! Momentum Neutral Transfert Coefficients (should be a constant) 1154 1341 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 1155 1342 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7 1156 zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details)1343 zCdn_ice = zCdn_skin_ice ! Eq. 7 1157 1344 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32) 1158 1345 1159 1346 ! Heat Neutral Transfert Coefficients 1160 zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) ) ! Eq. 50 + Eq. 52 (cf Lupkes email for details)1161 1347 zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) ) ! Eq. 50 + Eq. 52 1348 1162 1349 ! Atmospheric and Surface Variables 1163 1350 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1164 zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) )! saturation humidity over ocean [kg/kg]1165 zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )! saturation humidity over ice [kg/kg]1351 zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg] 1352 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg] 1166 1353 ! 1167 1354 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility … … 1169 1356 ! Virtual potential temperature [K] 1170 1357 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1171 zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice1358 zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1172 1359 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1173 1360 1174 1361 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1175 1362 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1176 1363 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1177 1364 1178 1365 ! Momentum and Heat Neutral Transfert Coefficients 1179 1366 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1180 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1181 1182 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead )1367 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1368 1369 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?) 1183 1370 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1184 z0i = z0_skin_ice ! over ice (cf Lupkes email for details)1371 z0i = z0_skin_ice ! over ice 1185 1372 IF( zrib_o <= 0._wp ) THEN 1186 1373 zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 10 … … 1191 1378 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1192 1379 ENDIF 1193 1380 1194 1381 IF( zrib_i <= 0._wp ) THEN 1195 1382 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 … … 1199 1386 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1200 1387 ENDIF 1201 1388 1202 1389 ! Momentum Transfert Coefficients (Eq. 38) 1203 Cd(ji,jj) = zCdn_skin_ice * zfmi + &1390 pcd(ji,jj) = zCdn_skin_ice * zfmi + & 1204 1391 & zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 1205 1392 1206 1393 ! Heat Transfert Coefficients (Eq. 49) 1207 Ch(ji,jj) = zChn_skin_ice * zfhi + &1394 pch(ji,jj) = zChn_skin_ice * zfhi + & 1208 1395 & zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 1209 1396 ! 1210 1397 END DO 1211 1398 END DO 1212 CALL lbc_lnk_multi( 'sbcblk', Cd, 'T', 1., Ch, 'T', 1. )1399 CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1., pch, 'T', 1. ) 1213 1400 ! 1214 1401 END SUBROUTINE Cdn10_Lupkes2015 -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r10069 r12154 1 1 MODULE sbcblk_algo_ecmwf 2 2 !!====================================================================== 3 !! *** MODULE sbcblk_algo_ecmwf *** 4 !! Computes turbulent components of surface fluxes 5 !! according to the method in IFS of the ECMWF model 6 !! 3 !! *** MODULE sbcblk_algo_ecmwf *** 4 !! Computes: 7 5 !! * bulk transfer coefficients C_D, C_E and C_H 8 6 !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed … … 10 8 !! => all these are used in bulk formulas in sbcblk.F90 11 9 !! 12 !! Using the bulk formulation/param. of IFS of ECMWF (cycle 31r2)10 !! Using the bulk formulation/param. of IFS of ECMWF (cycle 40r1) 13 11 !! based on IFS doc (avaible online on the ECMWF's website) 14 12 !! 13 !! Routine turb_ecmwf maintained and developed in AeroBulk 14 !! (https://github.com/brodeau/aerobulk) 15 15 !! 16 !! Routine turb_ecmwf maintained and developed in AeroBulk 17 !! (http://aerobulk.sourceforge.net/) 18 !! 19 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 16 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 20 17 !!---------------------------------------------------------------------- 21 18 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code … … 41 38 42 39 USE sbc_oce ! Surface boundary condition: ocean fields 40 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 41 USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB 43 42 44 43 IMPLICIT NONE 45 44 PRIVATE 46 45 47 PUBLIC :: TURB_ECMWF ! called by sbcblk.F9048 49 ! !! ECMWF own values for given constants, taken form IFS documentation...46 PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF 47 48 !! ECMWF own values for given constants, taken form IFS documentation... 50 49 REAL(wp), PARAMETER :: charn0 = 0.018 ! Charnock constant (pretty high value here !!! 51 50 ! ! => Usually 0.011 for moderate winds) 52 51 REAL(wp), PARAMETER :: zi0 = 1000. ! scale height of the atmospheric boundary layer...1 53 52 REAL(wp), PARAMETER :: Beta0 = 1. ! gustiness parameter ( = 1.25 in COAREv3) 54 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature...55 REAL(wp), PARAMETER :: Cp_dry = 1005.0 ! Specic heat of dry air, constant pressure [J/K/kg]56 REAL(wp), PARAMETER :: Cp_vap = 1860.0 ! Specic heat of water vapor, constant pressure [J/K/kg]57 53 REAL(wp), PARAMETER :: alpha_M = 0.11 ! For roughness length (smooth surface term) 58 54 REAL(wp), PARAMETER :: alpha_H = 0.40 ! (Chapter 3, p.34, IFS doc Cy31r1) 59 55 REAL(wp), PARAMETER :: alpha_Q = 0.62 ! 56 57 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 58 60 59 !!---------------------------------------------------------------------- 61 60 CONTAINS 62 61 63 SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu, & 64 & Cd, Ch, Ce , t_zu, q_zu, U_blk, & 65 & Cdn, Chn, Cen ) 66 !!---------------------------------------------------------------------------------- 67 !! *** ROUTINE turb_ecmwf *** 68 !! 69 !! 2015: L. Brodeau (brodeau@gmail.com) 70 !! 71 !! ** Purpose : Computes turbulent transfert coefficients of surface 72 !! fluxes according to IFS doc. (cycle 31) 73 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 74 !! 75 !! ** Method : Monin Obukhov Similarity Theory 62 63 SUBROUTINE sbcblk_algo_ecmwf_init(l_use_cs, l_use_wl) 64 !!--------------------------------------------------------------------- 65 !! *** FUNCTION sbcblk_algo_ecmwf_init *** 76 66 !! 77 67 !! INPUT : 78 68 !! ------- 69 !! * l_use_cs : use the cool-skin parameterization 70 !! * l_use_wl : use the warm-layer parameterization 71 !!--------------------------------------------------------------------- 72 LOGICAL , INTENT(in) :: l_use_cs ! use the cool-skin parameterization 73 LOGICAL , INTENT(in) :: l_use_wl ! use the warm-layer parameterization 74 INTEGER :: ierr 75 !!--------------------------------------------------------------------- 76 IF( l_use_wl ) THEN 77 ierr = 0 78 ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 79 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_wl & Hz_wl failed!' ) 80 dT_wl(:,:) = 0._wp 81 Hz_wl(:,:) = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars) 82 ENDIF 83 IF( l_use_cs ) THEN 84 ierr = 0 85 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 86 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' ) 87 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 88 ENDIF 89 END SUBROUTINE sbcblk_algo_ecmwf_init 90 91 92 93 SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 94 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 95 & Cdn, Chn, Cen, & 96 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 97 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 98 !!---------------------------------------------------------------------- 99 !! *** ROUTINE turb_ecmwf *** 100 !! 101 !! ** Purpose : Computes turbulent transfert coefficients of surface 102 !! fluxes according to IFS doc. (cycle 45r1) 103 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 104 !! Returns the effective bulk wind speed at zu to be used in the bulk formulas 105 !! 106 !! Applies the cool-skin warm-layer correction of the SST to T_s 107 !! if the net shortwave flux at the surface (Qsw), the downwelling longwave 108 !! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 109 !! are provided as (optional) arguments! 110 !! 111 !! INPUT : 112 !! ------- 113 !! * kt : current time step (starts at 1) 79 114 !! * zt : height for temperature and spec. hum. of air [m] 80 !! * zu : height for wind speed (generally 10m) [m] 81 !! * U_zu : scalar wind speed at 10m [m/s] 82 !! * sst : SST [K] 115 !! * zu : height for wind speed (usually 10m) [m] 83 116 !! * t_zt : potential air temperature at zt [K] 84 !! * ssq : specific humidity at saturation at SST [kg/kg]85 117 !! * q_zt : specific humidity of air at zt [kg/kg] 86 !! 118 !! * U_zu : scalar wind speed at zu [m/s] 119 !! * l_use_cs : use the cool-skin parameterization 120 !! * l_use_wl : use the warm-layer parameterization 121 !! 122 !! INPUT/OUTPUT: 123 !! ------------- 124 !! * T_s : always "bulk SST" as input [K] 125 !! -> unchanged "bulk SST" as output if CSWL not used [K] 126 !! -> skin temperature as output if CSWL used [K] 127 !! 128 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 129 !! -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True) 130 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False) 131 !! 132 !! OPTIONAL INPUT: 133 !! --------------- 134 !! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2] 135 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 136 !! * slp : sea-level pressure [Pa] 137 !! 138 !! OPTIONAL OUTPUT: 139 !! ---------------- 140 !! * pdT_cs : SST increment "dT" for cool-skin correction [K] 141 !! * pdT_wl : SST increment "dT" for warm-layer correction [K] 142 !! * pHz_wl : thickness of warm-layer [m] 87 143 !! 88 144 !! OUTPUT : … … 93 149 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 94 150 !! * q_zu : specific humidity of air // [kg/kg] 95 !! * U_blk : bulk wind at 10m [m/s] 96 !! 97 !! 98 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 99 !!---------------------------------------------------------------------------------- 151 !! * U_blk : bulk wind speed at zu [m/s] 152 !! 153 !! 154 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 155 !!---------------------------------------------------------------------------------- 156 INTEGER, INTENT(in ) :: kt ! current time step 100 157 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 101 158 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] 102 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst! sea surface temperature [Kelvin]159 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin] 103 160 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 104 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq! sea surface specific humidity [kg/kg]105 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity 161 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg] 162 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 163 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 164 LOGICAL , INTENT(in ) :: l_use_cs ! use the cool-skin parameterization 165 LOGICAL , INTENT(in ) :: l_use_wl ! use the warm-layer parameterization 107 166 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) 108 167 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens) … … 110 169 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 111 170 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg] 112 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind at 10m[m/s]171 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 113 172 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 173 ! 174 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 175 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] 176 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 177 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_cs 178 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K] 179 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 180 ! 115 181 INTEGER :: j_itt 116 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 117 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations 118 ! 119 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star, & 120 & dt_zu, dq_zu, & 121 & znu_a, & !: Nu_air, Viscosity of air 122 & Linv, & !: 1/L (inverse of Monin Obukhov length... 123 & z0, z0t, z0q 124 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 125 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 126 !!---------------------------------------------------------------------------------- 127 ! 128 ! Identical first gess as in COARE, with IFS parameter values though 129 ! 182 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 183 ! 184 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star 185 REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu 186 REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air 187 REAL(wp), DIMENSION(jpi,jpj) :: Linv !: 1/L (inverse of Monin Obukhov length... 188 REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q 189 ! 190 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST 191 ! 192 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 193 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 194 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90' 195 !!---------------------------------------------------------------------------------- 196 197 IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 198 130 199 l_zt_equal_zu = .FALSE. 131 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 132 133 200 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 201 202 !! Initializations for cool skin and warm layer: 203 IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 204 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 205 206 IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 207 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 208 209 IF( l_use_cs .OR. l_use_wl ) THEN 210 ALLOCATE ( zsst(jpi,jpj) ) 211 zsst = T_s ! backing up the bulk SST 212 IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction 213 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 214 ENDIF 215 216 217 ! Identical first gess as in COARE, with IFS parameter values though... 218 ! 134 219 !! First guess of temperature and humidity at height zu: 135 t_zu = MAX( t_zt , 0.0) ! who knows what's given on masked-continental regions...136 q_zu = MAX( q_zt , 1.e-6 ) ! "220 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... 221 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 137 222 138 223 !! Pot. temp. difference (and we don't want it to be 0!) 139 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.e-6), dt_zu )140 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.e-9), dq_zu )141 142 znu_a = visc_air(t_z t) ! Air viscosity (m^2/s) at zt given from temperature in (K)143 144 ztmp2 = 0.5 * 0.5! initial guess for wind gustiness contribution145 U_blk = SQRT(U_zu*U_zu + ztmp2) 146 147 ! z0 = 0.0001148 ztmp2 = 10000. ! optimization: ztmp2 == 1/z0149 ztmp0 = LOG(zu*ztmp2) 150 z tmp1 = LOG(10.*ztmp2)151 u_star = 0.035*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10)152 153 z0 = charn0*u_star*u_star/grav + 0.11*znu_a/u_star154 z0t = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ! WARNING: 1/z0t !224 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 225 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 226 227 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 228 229 U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 230 231 ztmp0 = LOG( zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 232 ztmp1 = LOG(10._wp*10000._wp) ! " " " 233 u_star = 0.035_wp*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 234 235 z0 = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 236 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 237 238 z0t = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 239 z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 155 240 156 241 Cd = (vkarmn/ztmp0)**2 ! first guess of Cd 157 242 158 ztmp0 = vkarmn*vkarmn/LOG(zt *z0t)/Cd159 160 ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk ) ! Ribu = Bulk Richardson number161 162 !! First estimate of zeta_u, depending on the stability, ie sign of Ribu(ztmp2):163 ztmp1 = 0.5 + SIGN( 0.5 , ztmp2 )243 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 244 245 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 246 247 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 248 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 164 249 func_m = ztmp0*ztmp2 ! temporary array !! 165 !! Ribu < 0 Ribu > 0 Beta = 1.25166 func_h = (1.-ztmp1)*(func_m/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! temporary array !!! func_h == zeta_u167 & + ztmp1*(func_m*(1. + 27./9.*ztmp2/ztmp0))250 func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! BRN < 0 ! temporary array !!! func_h == zeta_u 251 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! BRN > 0 252 !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 168 253 169 254 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 170 ztmp0 = vkarmn/(LOG(zu*z0t) - psi_h_ecmwf(func_h))171 172 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h))255 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 256 257 u_star = MAX ( U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on) 173 258 t_star = dt_zu*ztmp0 174 259 q_star = dq_zu*ztmp0 175 260 176 ! What 's needto be done if zt /= zu:261 ! What needs to be done if zt /= zu: 177 262 IF( .NOT. l_zt_equal_zu ) THEN 178 !179 263 !! First update of values at zu (or zt for wind) 180 264 ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu) ! zt*func_h/zu == zeta_t 181 ztmp1 = log(zt/zu) + ztmp0265 ztmp1 = LOG(zt/zu) + ztmp0 182 266 t_zu = t_zt - t_star/vkarmn*ztmp1 183 267 q_zu = q_zt - q_star/vkarmn*ztmp1 184 q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 185 186 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 187 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 268 q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 188 269 ! 270 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 271 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 189 272 ENDIF 190 273 … … 194 277 195 278 !! First guess of inverse of Monin-Obukov length (1/L) : 196 ztmp0 = (1. + rctv0*q_zu) ! the factor to apply to temp. to get virt. temp... 197 Linv = grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / ( u_star*u_star * t_zu*ztmp0 ) 279 Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 198 280 199 281 !! Functions such as u* = U_blk*vkarmn/func_m 200 ztmp1 = zu + z0 201 ztmp0 = ztmp1*Linv 202 func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv) 203 func_h = LOG(ztmp1*z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(1./z0t*Linv) 204 282 ztmp0 = zu*Linv 283 func_m = LOG(zu) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 284 func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 205 285 206 286 !! ITERATION BLOCK 207 !! ***************208 209 287 DO j_itt = 1, nb_itt 210 288 211 289 !! Bulk Richardson Number at z=zu (Eq. 3.25) 212 ztmp0 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk)290 ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 213 291 214 292 !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 215 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3, p.33, IFS doc - Cy31r1 293 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 294 !! Note: it is slightly different that the L we would get with the usual 295 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) 216 296 217 297 !! Update func_m with new Linv: 218 ztmp1 = zu + z0 219 func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp1*Linv) + psi_m_ecmwf(z0*Linv) 298 func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zu*Linv) + psi_m_ecmwf(z0*Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu! 220 299 221 300 !! Need to update roughness lengthes: … … 223 302 ztmp2 = u_star*u_star 224 303 ztmp1 = znu_a/u_star 225 z0 = alpha_M*ztmp1 + charn0*ztmp2/grav 226 z0t = alpha_H*ztmp1 ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 227 z0q = alpha_Q*ztmp1 228 229 !! Update wind at 10m taking into acount convection-related wind gustiness: 230 ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!! 231 ztmp2 = ztmp2 * (MAX(-zi0*Linv/vkarmn,0.))**(2./3.) ! => w*^2 (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 232 !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..) 233 U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2) ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1 304 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 305 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 306 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) 307 308 !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8) 309 ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 310 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 311 U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 234 312 ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 235 313 … … 238 316 !! as well the air-sea differences: 239 317 IF( .NOT. l_zt_equal_zu ) THEN 240 241 318 !! Arrays func_m and func_h are free for a while so using them as temporary arrays... 242 func_h = psi_h_ecmwf( (zu+z0)*Linv) ! temporary array !!!243 func_m = psi_h_ecmwf( (zt+z0)*Linv) ! temporary array !!!319 func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!! 320 func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!! 244 321 245 322 ztmp2 = psi_h_ecmwf(z0t*Linv) 246 323 ztmp0 = func_h - ztmp2 247 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0t) - ztmp0)324 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 248 325 t_star = dt_zu*ztmp1 249 326 ztmp2 = ztmp0 - func_m + ztmp2 … … 253 330 ztmp2 = psi_h_ecmwf(z0q*Linv) 254 331 ztmp0 = func_h - ztmp2 255 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0q) - ztmp0)332 ztmp1 = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0) 256 333 q_star = dq_zu*ztmp1 257 334 ztmp2 = ztmp0 - func_m + ztmp2 258 ztmp1 = log(zt/zu) + ztmp2335 ztmp1 = LOG(zt/zu) + ztmp2 259 336 q_zu = q_zt - q_star/vkarmn*ztmp1 260 261 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 262 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 263 264 END IF 337 ENDIF 265 338 266 339 !! Updating because of updated z0 and z0t and new Linv... 267 ztmp1 = zu + z0 268 ztmp0 = ztmp1*Linv 269 func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 270 func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 271 272 END DO 340 ztmp0 = zu*Linv 341 func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 342 func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 343 344 345 IF( l_use_cs ) THEN 346 !! Cool-skin contribution 347 348 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 349 & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0 350 351 CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst ) ! Qnsol -> ztmp1 352 353 T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 354 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 355 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 356 357 ENDIF 358 359 IF( l_use_wl ) THEN 360 !! Warm-layer contribution 361 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 362 & ztmp1, ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp2 363 CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 364 !! Updating T_s and q_s !!! 365 T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) ! 366 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 367 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 368 ENDIF 369 370 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 371 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 372 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 373 ENDIF 374 375 END DO !DO j_itt = 1, nb_itt 273 376 274 377 Cd = vkarmn*vkarmn/(func_m*func_m) 275 378 Ch = vkarmn*vkarmn/(func_m*func_h) 276 ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 277 Ce = vkarmn*vkarmn/(func_m*ztmp1) 278 279 ztmp1 = zu + z0 280 Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 281 Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 282 Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q)) 283 284 END SUBROUTINE TURB_ECMWF 379 ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 380 Ce = vkarmn*vkarmn/(func_m*ztmp2) 381 382 Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 383 Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 384 Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 385 386 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 387 IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 388 IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 389 390 IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 391 392 END SUBROUTINE turb_ecmwf 285 393 286 394 … … 294 402 !! and L is M-O length 295 403 !! 296 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)404 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 297 405 !!---------------------------------------------------------------------------------- 298 406 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf … … 302 410 REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 303 411 !!---------------------------------------------------------------------------------- 304 !305 412 DO jj = 1, jpj 306 413 DO ji = 1, jpi 307 414 ! 308 zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!):415 zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 309 416 ! 310 417 ! Unstable (Paulson 1970): 311 418 ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 312 zx = SQRT(ABS(1. - 16.*zzeta))313 ztmp = 1. + SQRT(zx)419 zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 420 ztmp = 1._wp + SQRT(zx) 314 421 ztmp = ztmp*ztmp 315 psi_unst = LOG( 0.125 *ztmp*(1.+ zx) ) &316 & -2. *ATAN( SQRT(zx) ) + 0.5*rpi422 psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) ) & 423 & -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 317 424 ! 318 425 ! Unstable: 319 426 ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 320 psi_stab = -2. /3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) &321 & - zzeta - 2. /3.*5./0.35427 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 428 & - zzeta - 2._wp/3._wp*5._wp/0.35_wp 322 429 ! 323 430 ! Combining: 324 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1325 ! 326 psi_m_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable327 & + stab * psi_stab ! (zzeta > 0) Stable431 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 432 ! 433 psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 434 & + stab * psi_stab ! (zzeta > 0) Stable 328 435 ! 329 436 END DO 330 437 END DO 331 !332 438 END FUNCTION psi_m_ecmwf 333 439 334 440 335 441 FUNCTION psi_h_ecmwf( pzeta ) 336 442 !!---------------------------------------------------------------------------------- … … 342 448 !! and L is M-O length 343 449 !! 344 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)450 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 345 451 !!---------------------------------------------------------------------------------- 346 452 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf … … 354 460 DO ji = 1, jpi 355 461 ! 356 zzeta = MIN(pzeta(ji,jj) , 5. ) ! Very stable conditions (L positif and big!):357 ! 358 zx = ABS(1. - 16.*zzeta)**.25 ! this is actually (1/phi_m)**2 !!!462 zzeta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!): 463 ! 464 zx = ABS(1._wp - 16._wp*zzeta)**.25 ! this is actually (1/phi_m)**2 !!! 359 465 ! ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 360 466 ! Unstable (Paulson 1970) : 361 psi_unst = 2. *LOG(0.5*(1.+ zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1467 psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 362 468 ! 363 469 ! Stable: 364 psi_stab = -2. /3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1365 & - ABS(1. + 2./3.*zzeta)**1.5 - 2./3.*5./0.35 + 1.470 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 471 & - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 366 472 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 367 473 ! 368 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1369 ! 370 ! 371 psi_h_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable372 & + stab * psi_stab ! (zzeta > 0) Stable474 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 475 ! 476 ! 477 psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 478 & + stab * psi_stab ! (zzeta > 0) Stable 373 479 ! 374 480 END DO 375 481 END DO 376 !377 482 END FUNCTION psi_h_ecmwf 378 483 379 380 FUNCTION Ri_bulk( pz, ptz, pdt, pqz, pdq, pub )381 !!----------------------------------------------------------------------------------382 !! Bulk Richardson number (Eq. 3.25 IFS doc)383 !!384 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)385 !!----------------------------------------------------------------------------------386 REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk !387 <