diff --git a/extern/fabm b/extern/fabm index 2e580282..717e63e3 160000 --- a/extern/fabm +++ b/extern/fabm @@ -1 +1 @@ -Subproject commit 2e5802828fe242d3144740332224db424e142cde +Subproject commit 717e63e3f4922e0975d4eebb20af8998ee194d43 diff --git a/extern/flexout b/extern/flexout index f129ad21..4c141f25 160000 --- a/extern/flexout +++ b/extern/flexout @@ -1 +1 @@ -Subproject commit f129ad2181576624ae929331388d0a69bc574f7f +Subproject commit 4c141f253dbc1e37657fbcf00be0c2da611f9ec0 diff --git a/extern/stim b/extern/stim index a8d2346f..bd929672 160000 --- a/extern/stim +++ b/extern/stim @@ -1 +1 @@ -Subproject commit a8d2346f08b61012d4ba015eb5a9d3d1a86305df +Subproject commit bd9296724cd2911eb30fc2ef6c9d17dd70b76a5c diff --git a/src/airsea/airsea.F90 b/src/airsea/airsea.F90 index 34b0e832..af4f8d75 100644 --- a/src/airsea/airsea.F90 +++ b/src/airsea/airsea.F90 @@ -65,6 +65,7 @@ module airsea_driver ! ! Meteorological forcing variables integer, public :: hum_method + integer, public :: longwave_type character(len=PATH_MAX) :: meteo_file type (type_scalar_input), public, target :: u10,v10 type (type_scalar_input), public, target :: airp,airt @@ -76,10 +77,10 @@ module airsea_driver ! ! surface shortwave radiation ! and surface heat flux (W/m^2) - type (type_scalar_input), public, target :: I_0, ql + type (type_scalar_input), public, target :: I_0 REALTYPE, public, target :: albedo - type (type_scalar_input), public, target :: heat - REALTYPE, public :: qe,qh + type (type_scalar_input), public, target :: heat, ql_ + REALTYPE, public :: qe,qh,ql ! surface stress components (Pa) REALTYPE, public, target :: tx,ty @@ -341,7 +342,7 @@ subroutine init_airsea_nml(namlst, fn) call tx_%configure(method=momentum_method, path=momentumflux_file, index=1, constant_value=const_tx) call ty_%configure(method=momentum_method, path=momentumflux_file, index=2, constant_value=const_ty) call heat%configure(method=heat_method, path=heatflux_file, index=1, scale_factor=shf_factor, constant_value=const_heat) - call ql%configure(method=back_radiation_method, path=back_radiation_file, index=1) + call ql_%configure(method=back_radiation_method, path=back_radiation_file, index=1) call sst_obs%configure(method=sst_method, path=sst_file, index=1) call sss%configure(method=sss_method, path=sss_file, index=1) call precip%configure(method=precip_method, path=precip_file, index=1, scale_factor=precip_factor, constant_value=const_precip) @@ -491,12 +492,24 @@ subroutine init_airsea_yaml() call branch%get(calc_evaporation, 'calc_evaporation', 'calculate evaporation from meteorological conditions', & default=.false.) + !jpnote: before cherry pick call branch%get(I_0, 'swr', 'shortwave radiation', 'W/m^2', & minimum=0._rk,default=0._rk, extra_options=(/option(3, 'from time, location and cloud cover', 'calculate')/)) - call branch%get(ql, 'longwave_radiation', 'net longwave radiation', 'W/m^2', & + !call branch%get(ql, 'longwave_radiation', 'net longwave radiation', 'W/m^2', & + ! default=0._rk, method_file=0, method_constant=method_unsupported, & + !extra_options=(/option(CLARK, 'Clark et al. (1974)', 'Clark'), option(HASTENRATH_LAMB, 'Hastenrath and Lamb (1978)', 'Hastenrath_Lamb'), option(BIGNAMI, 'Bignami et al. (1995)', 'Bignami'), option(BERLIAND_BERLIAND, 'Berliand and Berliand (1952)', 'Berliand_Berliand'), option(JOSEY1, 'Josey et al. (2003) - 1', 'Josey1'), option(JOSEY2, 'Josey et al. (2003) - 2', 'Josey2')/), default_method=CLARK) +!jpnote: after cherry pick 1 + !call branch%get(ql, 'longwave_radiation', 'longwave radiation', 'W/m^2', & + !call twig%get(ssuv_method, 'ssuv_method', 'wind treatment', & + !options=(/option(0, 'use absolute wind speed'), option(1, 'use wind speed relative to current velocity')/), default=1, display=display_advanced) +!jpnote: after cherry pick 2 + call branch%get(ql_, 'longwave_radiation', 'longwave radiation', 'W/m^2', & default=0._rk, method_file=0, method_constant=method_unsupported, & - extra_options=(/option(CLARK, 'Clark et al. (1974)', 'Clark'), option(HASTENRATH_LAMB, 'Hastenrath and Lamb (1978)', 'Hastenrath_Lamb'), option(BIGNAMI, 'Bignami et al. (1995)', 'Bignami'), option(BERLIAND_BERLIAND, 'Berliand and Berliand (1952)', 'Berliand_Berliand'), option(JOSEY1, 'Josey et al. (2003) - 1', 'Josey1'), option(JOSEY2, 'Josey et al. (2003) - 2', 'Josey2')/), default_method=CLARK) + extra_options=(/option(1, 'Clark'), option(2, 'Hastenrath'), option(3, 'Bignami'), option(4, 'Berliand'), option(5, 'Josey-1'), option(6, 'Josey-2')/), default_method=1, pchild=leaf) + call leaf%get(longwave_type, 'type', 'longwave type from file', & + options=(/option(1, 'net longwave radiation'), option(2, 'incoming longwave radiation')/), default=1) + twig => branch%get_typed_child('albedo') call twig%get(albedo_method, 'method', 'method to compute albedo', & options=(/option(0, 'constant', 'constant'), option(PAYNE, 'Payne (1972)', 'Payne'), option(COGLEY, 'Cogley (1979)', 'Cogley')/), default=PAYNE) @@ -687,11 +700,11 @@ subroutine post_init_airsea(lat,lon) LEVEL4 'using Fairall et. all formulation' case default end select - LEVEL3 'net longwave radiation:' - select case (ql%method) + LEVEL3 'longwave radiation:' + select case (ql_%method) case(0) ! Read from file instead of calculating - call register_input(ql) - case(CLARK) + call register_input(ql_) !jpnote commmit 2 change + case(1) LEVEL4 'using Clark formulation' case(HASTENRATH_LAMB) LEVEL4 'using Hastenrath formulation' @@ -762,9 +775,11 @@ subroutine surface_fluxes_uvic(surface_temp,sensible,latent,longwave_radiation_v ! call humidity(hum_method,rh,airp,TTss-kelvin,airt) call humidity(hum_method,hum%value,airp%value,tw,airt%value) ! call longwave_radiation(longwave_radiation_method, & - ! lat,TTss,airt+kelvin,cloud,qb) - call longwave_radiation(ql%method, & - dlat,tw_k,ta_k,cloud%value,longwave_radiation_value) + ! lat,TTss,airt+kelvin,cloud,qb) + call longwave_radiation(ql_%method,longwave_type, & + dlat,tw_k,ta_k,cloud%value,ql_%value,longwave_radiation_value) + !call longwave_radiation(ql%method, & + ! dlat,tw_k,ta_k,cloud%value,longwave_radiation_value) !call airsea_fluxes(fluxes_method, & !TTss-kelvin,airt,u10,v10,precip,evap,tx,ty,qe,qh) call airsea_fluxes(fluxes_method, & @@ -792,10 +807,10 @@ subroutine surface_fluxes(surface_temp,sensible,latent,longwave_radiation) latent = qe #if 0 if (qe .lt. _ZERO_) then - STDERR 'Stefan# ',qh/qe +!KB STDERR 'Stefan# ',qh/qe end if #endif - longwave_radiation = ql%value + longwave_radiation = ql return end subroutine surface_fluxes !EOC @@ -1017,9 +1032,8 @@ subroutine flux_from_meteo(jul,secs) cloud1 = cloud2 call humidity(hum_method,hum,airp,tw,ta) - if (ql%method .gt. 0) then - call longwave_radiation(ql%method, & - dlat,tw_k,ta_k,cloud,ql) + call longwave_radiation(ql_%method,longwave_type, & + dlat,tw_k,ta_k,cloud,ql_%value,ql) end if #if 0 call airsea_fluxes(fluxes_method,rain_impact,calc_evaporation, & @@ -1028,7 +1042,7 @@ subroutine flux_from_meteo(jul,secs) call airsea_fluxes(fluxes_method, & tw,ta,u10%value-ssu,v10%value-ssv,precip%value,evap,tx2,ty2,qe,qh) #endif - h2 = ql%value+qe+qh + h2 = ql+qe+qh cloud2 = cloud%value if (init_saved_vars) then @@ -1069,13 +1083,11 @@ subroutine flux_from_meteo(jul,secs) end if call humidity(hum_method,hum%value,airp%value,tw,ta) - if (ql%method .gt. 0) then - call longwave_radiation(ql%method, & - dlat,tw_k,ta_k,cloud%value,ql%value) - endif + call longwave_radiation(ql_%method,longwave_type, & + dlat,tw_k,ta_k,cloud%value,ql_%value,ql) call airsea_fluxes(fluxes_method, & tw,ta,u10%value-ssu,v10%value-ssv,precip%value,evap,tx_%value,ty_%value,qe,qh) - heat%value = (ql%value+qe+qh) + heat%value = (ql+qe+qh) #endif w = sqrt((u10%value-ssu)*(u10%value-ssu)+(v10%value-ssv)*(v10%value-ssv)) diff --git a/src/airsea/airsea_variables.F90 b/src/airsea/airsea_variables.F90 index 813273a7..7773fa70 100644 --- a/src/airsea/airsea_variables.F90 +++ b/src/airsea/airsea_variables.F90 @@ -23,7 +23,8 @@ module airsea_variables REALTYPE, public, parameter :: emiss=0.97 REALTYPE, public, parameter :: bolz=5.67e-8 REALTYPE, public, parameter :: kelvin=273.15 - REALTYPE, public, parameter :: const06=0.62198 + REALTYPE, public, parameter :: const06=0.62198 ! molecular weight ratio between water and dry air + ! 18.01528 g/mol H2O, 28.97 g/mol dry air REALTYPE, public, parameter :: rgas = 287.1 ! REALTYPE, public, parameter :: g = 9.81 ! [m/s2] REALTYPE, public, parameter :: rho_0 = 1025. ! [kg/m3] @@ -46,6 +47,7 @@ module airsea_variables integer, public, parameter :: COGLEY=2 ! Longwave radiation + integer, public, parameter :: FILE = 0 ! From file integer, public, parameter :: CLARK = 1 ! Clark et al, 1974 integer, public, parameter :: HASTENRATH_LAMB = 2 ! Hastenrath and Lamb, 1978 integer, public, parameter :: BIGNAMI = 3 ! Bignami et al., 1995 - Medsea diff --git a/src/airsea/longwave_radiation.F90 b/src/airsea/longwave_radiation.F90 index cbf32712..c98e086f 100644 --- a/src/airsea/longwave_radiation.F90 +++ b/src/airsea/longwave_radiation.F90 @@ -5,13 +5,14 @@ ! !ROUTINE: Calculate the net longwave radiation \label{sec:back-rad} ! ! !INTERFACE: - subroutine longwave_radiation(method,dlat,tw,ta,cloud,ql) + subroutine longwave_radiation(method,type,dlat,tw,ta,cloud,qlobs,ql) ! ! !DESCRIPTION: ! ! Here, the net longwave radiation is calculated by means of one out ! of six methods, which depend on the value given to the parameter ! {\tt method}: +! {\tt method}=0: read observations from a file ! {\tt method}=1: \cite{Clarketal74}, ! {\tt method}=2: \cite{HastenrathLamb78}, ! {\tt method}=3: \cite{Bignamietal95}, @@ -23,20 +24,29 @@ subroutine longwave_radiation(method,dlat,tw,ta,cloud,ql) ! !USES: use airsea_variables, only: emiss,bolz use airsea_variables, only: ea,qa - use airsea_variables, only: CLARK, HASTENRATH_LAMB, BIGNAMI, BERLIAND_BERLIAND, JOSEY1, JOSEY2 + use airsea_variables, only: FILE, CLARK, HASTENRATH_LAMB, BIGNAMI, BERLIAND_BERLIAND, JOSEY1, JOSEY2 IMPLICIT NONE ! ! !INPUT PARAMETERS: - integer, intent(in) :: method - REALTYPE, intent(in) :: dlat,tw,ta,cloud + integer, intent(in) :: method,type + REALTYPE, intent(in) :: dlat,tw,ta,cloud,qlobs ! ! !OUTPUT PARAMETERS: - REALTYPE, intent(out) :: ql + REALTYPE, intent(inout) :: ql ! ! !REVISION HISTORY: ! Original author(s): Adolf Stips, Hans Burchard & Karsten Bolding ! ! !LOCAL VARIABLES: +!cherry pick -- > jpnote to get rid of and add to airseavariables + !integer, parameter :: from_file=0 + !integer, parameter :: clark=1 ! Clark et al, 1974 + !integer, parameter :: hastenrath=2 ! Hastenrath and Lamb, 1978 + !integer, parameter :: bignami=3 ! Bignami et al., 1995 - Medsea + !integer, parameter :: berliand=4 ! Berliand and Berliand, 1952 - ROMS + !integer, parameter :: josey1=5 ! Josey 2003, (J1,9) + !integer, parameter :: josey2=6 ! Josey 2003, (J2,14) + REALTYPE, parameter, dimension(91) :: cloud_correction_factor = (/ & 0.497202, 0.501885, 0.506568, 0.511250, 0.515933, & 0.520616, 0.525299, 0.529982, 0.534665, 0.539348, & @@ -60,7 +70,6 @@ subroutine longwave_radiation(method,dlat,tw,ta,cloud,ql) REALTYPE :: ccf REALTYPE :: x1,x2,x3 -! !EOP !----------------------------------------------------------------------- !BOC @@ -68,6 +77,16 @@ subroutine longwave_radiation(method,dlat,tw,ta,cloud,ql) ccf= cloud_correction_factor(nint(abs(dlat))+1) select case(method) + !before + !case(CLARK) + !after + case(FILE) + select case(type) + case(1) + ql=qlobs + case(2) + ql = qlobs-bolz*emiss*(tw**4) + end select case(CLARK) ! Clark et al. (1974) formula. ! unit of ea is Pascal, must hPa diff --git a/src/fabm/gotm_fabm.F90 b/src/fabm/gotm_fabm.F90 index 2cb3a273..8945afa3 100644 --- a/src/fabm/gotm_fabm.F90 +++ b/src/fabm/gotm_fabm.F90 @@ -105,6 +105,8 @@ module gotm_fabm bioshade_feedback,bioalbedo_feedback,biodrag_feedback, & freshwater_impact,salinity_relaxation_to_freshwater_flux, & save_inputs, no_surface + + ! Arrays for work, vertical movement, and cross-boundary fluxes REALTYPE,allocatable,dimension(:,:) :: ws @@ -250,7 +252,7 @@ subroutine configure_gotm_fabm(cfg) ! ! !INPUT PARAMETERS: class (type_settings), intent(inout) :: cfg - type (type_settings), pointer :: branch + type (type_settings), pointer :: branch, twig ! ! !REVISION HISTORY: ! Original author(s): Jorn Bruggeman @@ -266,6 +268,7 @@ subroutine configure_gotm_fabm(cfg) ! Initialize all namelist variables to reasonable default values. call cfg%get(fabm_calc, 'use', 'enable FABM', & default=.false.) + call cfg%get(freshwater_impact, 'freshwater_impact', 'enable dilution/concentration by precipitation/evaporation', & default=.true.) ! disable to check mass conservation branch => cfg%get_child('feedbacks', 'feedbacks to physics') @@ -277,6 +280,7 @@ subroutine configure_gotm_fabm(cfg) default=.false.) call cfg%get(repair_state, 'repair_state', 'clip state to minimum/maximum boundaries', & default=.false.) + branch => cfg%get_child('numerics', display=display_advanced) call branch%get(ode_method, 'ode_method', 'time integration scheme applied to source terms', & options=(/ option(1, 'Forward Euler', 'FE'), option(2, 'Runge-Kutta 2', 'RK2'), option(3, 'Runge-Kutta 4', 'RK4'), & @@ -306,6 +310,8 @@ subroutine configure_gotm_fabm(cfg) options=(/option(-1, 'auto-detect (prefer fabm.yaml)', 'auto'), option(0, 'fabm.nml', 'nml'), option(1, 'fabm.yaml', 'yaml')/), & default=-1, display=display_advanced) + + LEVEL2 'done.' end subroutine configure_gotm_fabm @@ -699,6 +705,8 @@ subroutine init_var_gotm_fabm(nlev) ! and link it to FABM. decimal_yearday = _ZERO_ call model%link_scalar(standard_variables%number_of_days_since_start_of_the_year,decimal_yearday) + call model%link_scalar(standard_variables%timestep,dt) !jpnote added + allocate(Qsour(0:nlev),stat=rc) if (rc /= 0) stop 'init_var_gotm_fabm(): Error allocating (Qsour)' diff --git a/src/gotm/gotm.F90 b/src/gotm/gotm.F90 index 7ddc9724..d3542a7e 100644 --- a/src/gotm/gotm.F90 +++ b/src/gotm/gotm.F90 @@ -65,7 +65,6 @@ module gotm use stim_variables, only: ice_hs,ice_hi,ice_uvic_topmelt,ice_uvic_topgrowth,ice_uvic_termelt,ice_uvic_botmelt,ice_uvic_botgrowth,ice_uvic_tb,ice_uvic_ts use stim_variables, only: ice_uvic_Fs,ice_uvic_Ff,ice_uvic_parb,ice_uvic_parui,ice_uvic_Amelt - use turbulence, only: turb_method use turbulence, only: init_turbulence,post_init_turbulence,do_turbulence use turbulence, only: num,nuh,nus @@ -569,6 +568,27 @@ subroutine initialize_gotm() call model_fabm%link_horizontal_data(standard_variables_fabm%bottom_depth,depth) call model_fabm%link_horizontal_data(standard_variables_fabm%bottom_depth_below_geoid,depth0) call model_fabm%link_horizontal_data(standard_variables_fabm%bottom_roughness_length,z0b) + + !ice vars--------- jpnote + + call model_fabm%link_horizontal_data(standard_variables_fabm%sea_ice_thickness,ice_hi) + call model_fabm%link_horizontal_data(standard_variables_fabm%snow_thickness,ice_hs) + call model_fabm%link_horizontal_data(standard_variables_fabm%topmelt,ice_uvic_topmelt) + call model_fabm%link_horizontal_data(standard_variables_fabm%f_melt,ice_uvic_Amelt) + call model_fabm%link_horizontal_data(standard_variables_fabm%topgrowth,ice_uvic_topgrowth) + call model_fabm%link_horizontal_data(standard_variables_fabm%termelt,ice_uvic_termelt) + call model_fabm%link_horizontal_data(standard_variables_fabm%tendency_of_sea_ice_thickness_due_to_thermodynamics_melt,ice_uvic_botmelt) + call model_fabm%link_horizontal_data(standard_variables_fabm%tendency_of_sea_ice_thickness_due_to_thermodynamics_grow,ice_uvic_botgrowth) + call model_fabm%link_horizontal_data(standard_variables_fabm%sea_ice_temperature,ice_uvic_tb) + call model_fabm%link_horizontal_data(standard_variables_fabm%surface_ice_temperature,ice_uvic_ts) + call model_fabm%link_horizontal_data(standard_variables_fabm%lowest_ice_layer_PAR,ice_uvic_parb) + call model_fabm%link_horizontal_data(standard_variables_fabm%under_ice_PAR,ice_uvic_parui) + + call model_fabm%link_interior_data(standard_variables_fabm%zonal_current,u(1:nlev)) !jpnote : changed from bulk data to interior data + call model_fabm%link_interior_data(standard_variables_fabm%meridional_current,v(1:nlev)) !jpnote: changed from bulk data to interior data + !------------------ + + if (fluxes_method /= 0) then call model_fabm%link_horizontal_data(standard_variables_fabm%surface_specific_humidity,qa) call model_fabm%link_horizontal_data(standard_variables_fabm%surface_air_pressure,airp%value) @@ -596,6 +616,9 @@ subroutine initialize_gotm() ! Initialize FABM initial state (this is done after the first call to do_input, ! to allow user-specified observed values to be used as initial state) if (fabm_calc) call init_gotm_fabm_state(nlev) + + + #endif if (restart) then @@ -800,7 +823,9 @@ subroutine integrate_gotm() ! update temperature and salinity if (sprof%method .ne. 0) then + call salinity(nlev,dt,cnpar,nus,gams,ice_uvic_Fs,ice_uvic_Ff) + endif if (tprof%method .ne. 0) then diff --git a/src/gotm/register_all_variables.F90 b/src/gotm/register_all_variables.F90 index 62711e54..44185e96 100644 --- a/src/gotm/register_all_variables.F90 +++ b/src/gotm/register_all_variables.F90 @@ -162,7 +162,8 @@ subroutine register_airsea_variables(nlev) call fm%register('I_0', 'W/m2', 'incoming short wave radiation', standard_name='', data0d=I_0%value, category='surface/heat_fluxes') call fm%register('qh', 'W/m2', 'sensible heat flux', standard_name='', data0d=qh, category='surface/heat_fluxes') call fm%register('qe', 'W/m2', 'latent heat flux', standard_name='', data0d=qe, category='surface/heat_fluxes') - call fm%register('ql', 'W/m2', 'net longwave radiation', standard_name='', data0d=ql%value, category='surface/heat_fluxes') + call fm%register('ql', 'W/m2', 'net longwave radiation', standard_name='', data0d=ql, category='surface/heat_fluxes') + call fm%register('qlobs', 'W/m2', 'longwave radiation (obs)', standard_name='', data0d=ql_%value, category='surface/heat_fluxes') call fm%register('heat', 'W/m2', 'net surface heat flux', standard_name='', data0d=heat%value, category='surface/heat_fluxes') call fm%register('tx', 'm2/s2', 'wind stress (x)', standard_name='', data0d=tx, category='surface') call fm%register('ty', 'm2/s2', 'wind stress (y)', standard_name='', data0d=ty, category='surface') diff --git a/src/meanflow/salinity.F90 b/src/meanflow/salinity.F90 index fd1fd4c6..49ce3d87 100644 --- a/src/meanflow/salinity.F90 +++ b/src/meanflow/salinity.F90 @@ -91,7 +91,7 @@ subroutine salinity(nlev,dt,cnpar,nus,gams,Fs,Ff) !jpnote passing in Fs and Ff f ! ! ! ice-ocean interaction salt and freshwater fluxes - REALTYPE, intent(in) :: Fs,Ff !jpnote added + REALTYPE, optional, intent(in) :: Fs,Ff ! ! !REVISION HISTORY: @@ -117,7 +117,9 @@ subroutine salinity(nlev,dt,cnpar,nus,gams,Fs,Ff) !jpnote passing in Fs and Ff f DiffBcup = Neumann DiffBcdw = Neumann DiffSup = -S(nlev)*(precip%value+evap+Ff)-Fs !NSnote, check signs -! DiffSup = -S(nlev)*(precip%value+evap) + + !DiffSup = -S(nlev)*(precip%value+evap) + DiffSdw = _ZERO_ AdvBcup = oneSided