diff --git a/.gitmodules b/.gitmodules index 44e38b7c..1e132023 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,14 +1,12 @@ [submodule "extern/fabm"] path = extern/fabm - url = https://github.com/fabm-model/fabm -# url = ../../fabm-model/fabm + url = https://github.com/PatrickFarnole/fabm [submodule "extern/flexout"] path = extern/flexout - url = https://github.com/BoldingBruggeman/flexout -# url = ../../BoldingBruggeman/flexout + url = https://github.com/PatrickFarnole/flexout [submodule "extern/stim"] path = extern/stim - url = https://github.com/BoldingBruggeman/stim + url = https://github.com/PatrickFarnole/stim [submodule "extern/CVMix-src"] path = extern/CVMix-src url = https://github.com/CVMix/CVMix-src.git diff --git a/extern/fabm b/extern/fabm index ab8a6ffc..717e63e3 160000 --- a/extern/fabm +++ b/extern/fabm @@ -1 +1 @@ -Subproject commit ab8a6ffc5ceecd6e467878fe7e3142adcc49ab3a +Subproject commit 717e63e3f4922e0975d4eebb20af8998ee194d43 diff --git a/extern/flexout b/extern/flexout index ff01d274..4c141f25 160000 --- a/extern/flexout +++ b/extern/flexout @@ -1 +1 @@ -Subproject commit ff01d2747564635eede8fe2943464b66d01b0b17 +Subproject commit 4c141f253dbc1e37657fbcf00be0c2da611f9ec0 diff --git a/extern/stim b/extern/stim index c24e322a..eb7d7af9 160000 --- a/extern/stim +++ b/extern/stim @@ -1 +1 @@ -Subproject commit c24e322a6de5c23f79da52f3e582757064cc2b94 +Subproject commit eb7d7af97dc77c4b124f2ca8de5e69bd8b4c835e diff --git a/src/airsea/airsea.F90 b/src/airsea/airsea.F90 index 45b30547..748d87aa 100644 --- a/src/airsea/airsea.F90 +++ b/src/airsea/airsea.F90 @@ -49,6 +49,7 @@ module airsea_driver public :: set_sst public :: set_ssuv public :: surface_fluxes + public :: surface_fluxes_uvic !jpnote added public :: integrated_fluxes #ifdef _PRINTSTATE_ public :: print_state_airsea @@ -64,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 @@ -75,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 @@ -262,7 +264,7 @@ subroutine init_airsea_nml(namlst, fn) REALTYPE :: const_tx,const_ty REALTYPE :: const_precip REALTYPE :: precip_factor - integer :: back_radiation_method + integer :: back_radiation_method integer :: heat_method namelist /airsea/ calc_fluxes, & @@ -340,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) @@ -354,6 +356,8 @@ subroutine init_airsea_nml(namlst, fn) stop 'init_airsea' 91 FATAL 'I could not read airsea namelist' stop 'init_airsea' +93 FATAL 'I could not open ',trim(meteo_file) + stop 'init_airsea' end subroutine init_airsea_nml !EOC @@ -492,10 +496,13 @@ subroutine init_airsea_yaml() 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', '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(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, 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) @@ -686,10 +693,10 @@ 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) + call register_input(ql_) case(CLARK) LEVEL4 'using Clark formulation' case(HASTENRATH_LAMB) @@ -735,6 +742,36 @@ subroutine post_init_airsea(lat,lon) end subroutine post_init_airsea !EOC +!----------------------------------------------------------------------- +!BOP + !jpnote added + subroutine surface_fluxes_uvic(surface_temp,sensible,latent,longwave_radiation_value) + ! + ! !USES: + IMPLICIT NONE + ! + ! !INPUT PARAMETERS: + REALTYPE, intent(in) :: surface_temp + ! !OUTPUT PARAMETERS: + REALTYPE, intent(out) :: sensible,latent,longwave_radiation_value + ! + ! LOCAL VARIABLES: + REALTYPE :: tw,tw_k,ta_k + !EOP + !----------------------------------------------------------------------- + !BOC + + tw = surface_temp-KELVIN + tw_k= surface_temp + ta_k = airt%value + KELVIN + + call humidity(hum_method,hum%value,airp%value,tw,airt%value) + call longwave_radiation(ql_%method,longwave_type, & + dlat,tw_k,ta_k,cloud%value,ql_%value,longwave_radiation_value) + call airsea_fluxes(fluxes_method, & + tw,airt%value,u10%value,v10%value,precip%value,evap,tx_%value,ty_%value,latent,sensible) + + end subroutine surface_fluxes_uvic !----------------------------------------------------------------------- !BOP subroutine surface_fluxes(surface_temp,sensible,latent,longwave_radiation) @@ -756,10 +793,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 @@ -981,10 +1018,9 @@ 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) - end if + call longwave_radiation(ql_%method,longwave_type, & + dlat,tw_k,ta_k,cloud,ql_%value,ql) + #if 0 call airsea_fluxes(fluxes_method,rain_impact,calc_evaporation, & tw,ta,u10%value-ssu,v10%value-ssv,precip%value,evap,tx2,ty2,qe,qh) @@ -992,7 +1028,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 @@ -1033,13 +1069,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..02e96e29 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,15 +24,15 @@ 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 @@ -60,7 +61,6 @@ subroutine longwave_radiation(method,dlat,tw,ta,cloud,ql) REALTYPE :: ccf REALTYPE :: x1,x2,x3 -! !EOP !----------------------------------------------------------------------- !BOC @@ -68,7 +68,14 @@ subroutine longwave_radiation(method,dlat,tw,ta,cloud,ql) ccf= cloud_correction_factor(nint(abs(dlat))+1) select case(method) - case(CLARK) + 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 ! Black body defect term, clouds, water vapor correction 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 5c99e9f7..5cea4065 100644 --- a/src/gotm/gotm.F90 +++ b/src/gotm/gotm.F90 @@ -49,18 +49,21 @@ module gotm use time use airsea_driver, only: init_airsea,post_init_airsea,do_airsea,clean_airsea - use airsea_driver, only: surface_fluxes + use airsea_driver, only: surface_fluxes, surface_fluxes_uvic !jpnote use airsea_driver, only: set_sst,set_ssuv,integrated_fluxes use airsea_driver, only: fluxes_method - use airsea_driver, only: wind=>w,tx,ty,I_0,cloud,heat,precip,evap,airp,albedo + use airsea_driver, only: wind=>w,tx,ty,hum,I_0,cloud,heat,precip,evap,airp,albedo use airsea_driver, only: bio_albedo,bio_drag_scale - use airsea_driver, only: u10, v10 + use airsea_driver, only: u10,v10,airt,sst,sss + use airsea_driver, only: ql,hum_method,fluxes_method use airsea_variables, only: qa,ta #ifdef _ICE_ use ice, only: init_ice, post_init_ice, do_ice, clean_ice, ice_cover - use stim_variables, only: Tice_surface,albedo_ice,transmissivity + use stim_variables, only: Tice_surface,albedo_ice,transmissivity,nilay,ice_uvic_Tice #endif + 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 @@ -412,6 +415,10 @@ subroutine initialize_gotm() call fm%register_dimension('lat',1,id=id_dim_lat) call fm%register_dimension('z',nlev,id=id_dim_z) call fm%register_dimension('zi',nlev+1,id=id_dim_zi) +#ifdef _ICE_ + call fm%register_dimension('zice',nilay+1,id=id_dim_zice) + call fm%register_dimension('dzice',nilay,id=id_dim_dzice) +#endif call fm%register_dimension('time',id=id_dim_time) call fm%initialize(prepend_by_default=(/id_dim_lon,id_dim_lat/),append_by_default=(/id_dim_time/)) @@ -561,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) @@ -588,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 @@ -684,6 +715,7 @@ subroutine integrate_gotm() !EOP ! ! !LOCAL VARIABLES: + integer(kind=timestepkind):: n integer :: i=0 @@ -743,7 +775,11 @@ subroutine integrate_gotm() #ifdef _ICE_ Qsw = I_0%value - call do_ice(h(nlev),dt,T(nlev),S(nlev),ta,precip%value,Qsw,surface_fluxes) + + call do_ice(h(nlev),dt,T(nlev),S(nlev),ta,precip%value,Qsw, & + surface_fluxes,surface_fluxes_uvic,julianday,secondsofday, & + I_0%value,airt%value,rho(nlev),rho_0,albedo,heat%value) + #endif ! reset some quantities @@ -781,7 +817,9 @@ subroutine integrate_gotm() ! update temperature and salinity if (sprof%method .ne. 0) then - call salinity(nlev,dt,cnpar,nus,gams) + + 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 00f0e50e..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') @@ -189,7 +190,17 @@ subroutine register_stim_variables(nlev) use stim_variables, only: Tice_surface,Tice,Tf use stim_variables, only: Hice, Hfrazil, dHis, dHib use stim_variables, only: surface_ice_energy,bottom_ice_energy + + use stim_variables, only: ice_hi,ice_hs,ice_uvic_hm,ice_uvic_ts,ice_uvic_tb,ice_uvic_swr_0,ice_uvic_precip_i,ice_uvic_sfall_i,ice_uvic_parb,ice_uvic_parui + use stim_variables, only: ice_uvic_Fh,ice_uvic_Fs,ice_uvic_Ff,ice_uvic_Sicebulk + use stim_variables, only: ice_uvic_topmelt,ice_uvic_botmelt,ice_uvic_termelt,ice_uvic_topgrowth,ice_uvic_botgrowth + use stim_variables, only: ice_uvic_Hmix,ice_uvic_Aice,ice_uvic_asnow,ice_uvic_Amelt,ice_uvic_dzi,ice_uvic_Cond,ice_uvic_rhoCp,ice_uvic_zi,ice_uvic_zice,ice_uvic_dzice,ice_uvic_Tice,ice_uvic_Sint,ice_uvic_Pari + use stim_variables, only: nilay use stim_variables, only: ocean_ice_flux + + + use stim_flato, only: simass,snmass,meltmass,qb + IMPLICIT NONE ! ! !INPUT PARAMETERS: @@ -206,6 +217,10 @@ subroutine register_stim_variables(nlev) call fm%register('T1', 'celsius', 'ice temperature (upper)', standard_name='', data0d=Tice(1), category='ice') call fm%register('T2', 'celsius', 'ice temperature (lower)', standard_name='', data0d=Tice(2), category='ice') end if + if (ice_model .eq. 4) then + call fm%register('T1', 'celsius', 'ice temperature (upper)', standard_name='', data0d=Tice(1), category='ice') + call fm%register('T2', 'celsius', 'ice temperature (lower)', standard_name='', data0d=Tice(2), category='ice') + end if call fm%register('Tf', 'celsius', 'ice freezing temperature', standard_name='', data0d=Tf, category='ice') call fm%register('Hice', 'm', 'ice thickness', standard_name='', data0d=Hice, category='ice') call fm%register('surface_ice_energy', 'J/m2', 'ice energy (surface)', standard_name='', data0d=surface_ice_energy, category='ice') @@ -215,6 +230,45 @@ subroutine register_stim_variables(nlev) call fm%register('dHis', 'm', 'ice growth (surface)', standard_name='', data0d=dHis, category='ice') call fm%register('dHib', 'm', 'ice growth (bottom)', standard_name='', data0d=dHib, category='ice') +! Flato +! Vars from ice_uvic.F90 + call fm%register('simass', 'kg/m2"', 'Ice mass per area', standard_name='', data0d=simass, category='ice') + call fm%register('snmass', 'kg/m2', 'Snow mass per area', standard_name='', data0d=snmass, category='ice') + call fm%register('meltmass', 'kg/m2', 'meltpond mass per area"', standard_name='', data0d=meltmass, category='ice') +! Public vars from ice.F90 + call fm%register('ice_hi', 'm', 'Ice thickness', standard_name='', data0d=ice_hi, category='ice') + call fm%register('ice_hs', 'm', 'Snow thickness', standard_name='', data0d=ice_hs, category='ice') + call fm%register('ice_hm', 'm', 'Meltpond thickness', standard_name='', data0d=ice_uvic_hm, category='ice') + call fm%register('ice_ts', 'K', 'ice/snow surface temperature', standard_name='', data0d=ice_uvic_ts, category='ice') + call fm%register('ice_tb', 'K', 'lowest ice layer temperature', standard_name='', data0d=ice_uvic_tb, category='ice') + call fm%register('swr_0', 'W/m^2', 'incidental shortwave radiation', standard_name='', data0d=ice_uvic_swr_0, category='ice') + call fm%register('precip_i', 'm/s', 'precipitation rate', standard_name='', data0d=ice_uvic_precip_i, category='ice') + call fm%register('sfall_i', 'm/s', 'snowfall rate', standard_name='', data0d=ice_uvic_sfall_i, category='ice') + call fm%register('ice_parb', 'W/m^2', 'lowest ice layer PAR', standard_name='', data0d=ice_uvic_parb, category='ice') + call fm%register('ice_parui', 'W/m^2', 'under ice PAR', standard_name='', data0d=ice_uvic_parui, category='ice') + call fm%register('Fh', 'W/m2', 'Heat flux at the interface', standard_name='', data0d=ice_uvic_fh, category='ice') + call fm%register('Fs', 'ppt m/s', 'Salt flux at the interface', standard_name='', data0d=ice_uvic_Fs, category='ice') + call fm%register('Ff', 'm/s', 'Freshwater flux at the interface', standard_name='', data0d=ice_uvic_Ff, category='ice') + call fm%register('Sice_bulk', 'ppt', 'Bulk sea ice salinity', standard_name='', data0d=ice_uvic_Sicebulk, category='ice') + call fm%register('topmelt', 'W/m2', 'Surface ice melt', standard_name='', data0d=ice_uvic_topmelt, category='ice') + call fm%register('botmelt', 'W/m2', 'Bottom ice melt', standard_name='', data0d=ice_uvic_botmelt, category='ice') + call fm%register('termelt', 'W/m2', 'Internal ice melt', standard_name='', data0d=ice_uvic_termelt, category='ice') + call fm%register('topgrowth', '', 'Surface ice growth (via submersion)', standard_name='', data0d=ice_uvic_topgrowth, category='ice') + call fm%register('botgrowth', '', 'Bottom ice growth', standard_name='', data0d=ice_uvic_botgrowth, category='ice') + call fm%register('Hmix', '', 'mixed layer heat storage', standard_name='', data0d=ice_uvic_Hmix, category='ice') + call fm%register('Aice', '', 'fraction of ice area that is bare ice', standard_name='', data0d=ice_uvic_Aice, category='ice') + call fm%register('Asnow', '', 'fraction of ice area that is snow covered', standard_name='', data0d=ice_uvic_Asnow, category='ice') + call fm%register('Amelt', '', 'fraction of ice area that is melt pond', standard_name='', data0d=ice_uvic_Amelt, category='ice') + call fm%register('dzi', 'm', 'Ice layer thickness', standard_name='', dimensions=(/id_dim_dzice/), data1d=ice_uvic_dzi(2:nilay+1), category='ice') + call fm%register('ice_uvic_zi', 'm', 'interface depth', standard_name='', dimensions=(/id_dim_zice/), data1d=ice_uvic_zi(2:nilay+2), category='ice') + call fm%register('zice', 'm', 'interface depth', standard_name='', dimensions=(/id_dim_zice/),data1d=ice_uvic_zice(1:nilay+1), category='ice') + call fm%register('dzice', 'm', 'interface depth', standard_name='', dimensions=(/id_dim_dzice/), data1d=ice_uvic_dzice(1:nilay), category='ice') + call fm%register('Cond', 'W m-1 kg-1', 'Ice conductivity', standard_name='',dimensions=(/id_dim_dzice/),data1d=ice_uvic_Cond(2:nilay+1), category='ice') + call fm%register('rhoCp', '10^6 J m-3 K-1', 'volum heat capacity', standard_name='',dimensions=(/id_dim_dzice/), data1d=ice_uvic_rhoCp(2:nilay+1), category='ice') + call fm%register('Tice', 'Celsius', 'Ice slab temperature', standard_name='',dimensions=(/id_dim_zice/), data1d=ice_uvic_Tice(2:nilay+2), category='ice') + call fm%register('Sint', 'W m-3', 'Ice SW internal heat', standard_name='',dimensions=(/id_dim_zice/),data1d=ice_uvic_Sint(2:nilay+2), category='ice') + call fm%register('Pari', 'W m-2', 'PAR radiation in ice', standard_name='',dimensions=(/id_dim_zice/), data1d=ice_uvic_Pari(2:nilay+2), category='ice') + return end subroutine register_stim_variables !EOC diff --git a/src/meanflow/salinity.F90 b/src/meanflow/salinity.F90 index 09a88082..49ce3d87 100644 --- a/src/meanflow/salinity.F90 +++ b/src/meanflow/salinity.F90 @@ -5,7 +5,7 @@ ! !ROUTINE: The salinity equation \label{sec:salinity} ! ! !INTERFACE: - subroutine salinity(nlev,dt,cnpar,nus,gams) + subroutine salinity(nlev,dt,cnpar,nus,gams,Fs,Ff) !jpnote passing in Fs and Ff from ice ! ! !DESCRIPTION: ! This subroutine computes the balance of salinity in the form @@ -89,6 +89,11 @@ subroutine salinity(nlev,dt,cnpar,nus,gams) ! non-local salinity flux (psu m/s) REALTYPE, intent(in) :: gams(0:nlev) ! +! +! ice-ocean interaction salt and freshwater fluxes + REALTYPE, optional, intent(in) :: Fs,Ff + +! ! !REVISION HISTORY: ! Original author(s): Hans Burchard & Karsten Bolding ! @@ -111,7 +116,10 @@ subroutine salinity(nlev,dt,cnpar,nus,gams) ! set boundary conditions DiffBcup = Neumann DiffBcdw = Neumann - DiffSup = -S(nlev)*(precip%value+evap) + DiffSup = -S(nlev)*(precip%value+evap+Ff)-Fs !NSnote, check signs + + !DiffSup = -S(nlev)*(precip%value+evap) + DiffSdw = _ZERO_ AdvBcup = oneSided