diff --git a/.github/workflows/gotm.yml b/.github/workflows/gotm.yml index 9bb033cd..337a99ec 100644 --- a/.github/workflows/gotm.yml +++ b/.github/workflows/gotm.yml @@ -15,7 +15,9 @@ jobs: runs-on: ${{ matrix.os }} steps: - name: Install compiler - run: sudo apt-get install gfortran-${{ matrix.version }} + run: | + sudo apt-get update + sudo apt-get install gfortran-${{ matrix.version }} - name: Clone gotm uses: actions/checkout@v2 with: diff --git a/extern/fabm b/extern/fabm index ab8a6ffc..c069e29d 160000 --- a/extern/fabm +++ b/extern/fabm @@ -1 +1 @@ -Subproject commit ab8a6ffc5ceecd6e467878fe7e3142adcc49ab3a +Subproject commit c069e29db82da530d9e8abf1bd2d442572e67da4 diff --git a/extern/stim b/extern/stim index fe0735ce..990f67d2 160000 --- a/extern/stim +++ b/extern/stim @@ -1 +1 @@ -Subproject commit fe0735ce55d96eafb9bc2d8b4063cc804dc72cde +Subproject commit 990f67d244af98094d6a76ad8956df8bc5cfbd90 diff --git a/src/fabm/gotm_fabm.F90 b/src/fabm/gotm_fabm.F90 index 40bc5f87..8945afa3 100644 --- a/src/fabm/gotm_fabm.F90 +++ b/src/fabm/gotm_fabm.F90 @@ -91,7 +91,7 @@ module gotm_fabm module procedure register_scalar_observation end interface - type (type_bulk_variable_id), save :: temp_id,salt_id,rho_id,h_id,swr_id,par_id,pres_id + type (type_bulk_variable_id), save :: temp_id,salt_id,rho_id,h_id,swr_id,par_id,pres_id,nuh_id type (type_horizontal_variable_id),save :: lon_id,lat_id,windspeed_id,par_sf_id,cloud_id,taub_id,swr_sf_id ! Variables to hold time spent on advection, diffusion, sink/source terms. @@ -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 @@ -113,7 +115,7 @@ module gotm_fabm integer,allocatable, dimension(:) :: posconc ! Arrays for environmental variables not supplied externally. - REALTYPE,allocatable,dimension(:),target :: par,pres,swr,k_par,z + REALTYPE,allocatable,dimension(:),target :: par,pres,swr,k_par,z,nuh_ct ! External variables REALTYPE :: dt,dt_eff ! External and internal time steps @@ -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 @@ -477,6 +483,7 @@ subroutine init_gotm_fabm(nlev,dt,field_manager) par_id = model%get_bulk_variable_id(standard_variables%downwelling_photosynthetic_radiative_flux) swr_id = model%get_bulk_variable_id(standard_variables%downwelling_shortwave_flux) pres_id = model%get_bulk_variable_id(standard_variables%pressure) + nuh_id = model%get_bulk_variable_id(type_interior_standard_variable(name='vertical_tracer_diffusivity', units='m s-1')) lon_id = model%get_horizontal_variable_id(standard_variables%longitude) lat_id = model%get_horizontal_variable_id(standard_variables%latitude) windspeed_id = model%get_horizontal_variable_id(standard_variables%wind_speed) @@ -678,6 +685,15 @@ subroutine init_var_gotm_fabm(nlev) call model%link_interior_data(pres_id,pres) end if + ! Allocate array for turbulent diffusivity at centers if necessary. + ! This will be calculated from diffusivity at interfaces during each time step. + if (model%variable_needs_values(nuh_id)) then + allocate(nuh_ct(1:nlev),stat=rc) + if (rc /= 0) stop 'init_var_gotm_fabm(): Error allocating (nuh_ct)' + nuh_ct = _ZERO_ + call model%link_interior_data(nuh_id, nuh_ct) + end if + ! Allocate array for local depth (below water surface). ! This will be calculated from layer depths. allocate(z(1:nlev),stat=rc) @@ -689,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)' @@ -1144,6 +1162,12 @@ subroutine calculate_derived_input(nlev) pres(1:nlev) = p0 + pres(1:nlev)*gravity/10000 end if + if (allocated(nuh_ct)) then + do i=1,nlev + nuh_ct(i) = 0.5_gotmrk * (curnuh(i - 1) + curnuh(i)) + end do + end if + ! Calculate local depth below surface from layer height ! Used internally to compute light field, and may be used by ! biogeochemical models as well. @@ -1669,11 +1693,12 @@ subroutine clean_gotm_fabm if (allocated(change_in_total)) deallocate(change_in_total) ! Deallocate arrays with internally computed environmental variables. - if (allocated(par)) deallocate(par) - if (allocated(k_par)) deallocate(k_par) - if (allocated(swr)) deallocate(swr) - if (allocated(pres)) deallocate(pres) - if (allocated(z)) deallocate(z) + if (allocated(par)) deallocate(par) + if (allocated(k_par)) deallocate(k_par) + if (allocated(swr)) deallocate(swr) + if (allocated(pres)) deallocate(pres) + if (allocated(z)) deallocate(z) + if (allocated(nuh_ct)) deallocate(nuh_ct) ! Reset all module-level pointers nullify(nuh) diff --git a/src/gotm/gotm.F90 b/src/gotm/gotm.F90 index 26136dac..9ae42abc 100644 --- a/src/gotm/gotm.F90 +++ b/src/gotm/gotm.F90 @@ -61,9 +61,9 @@ module gotm #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,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 -#endif use turbulence, only: turb_method use turbulence, only: init_turbulence,post_init_turbulence,do_turbulence @@ -415,11 +415,9 @@ 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) - ! call fm%register_dimension('z1',nilay,id=id_dim_z1) !jpnote commented #ifdef _ICE_ -!when compiling with stim off we dont have access to stim variables - call fm%register_dimension('zice',nilay+1,id=id_dim_zice) !jpnote - call fm%register_dimension('dzice',nilay,id=id_dim_dzice) !jpnote + 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/)) @@ -570,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) @@ -597,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 @@ -694,8 +716,6 @@ subroutine integrate_gotm() ! ! !LOCAL VARIABLES: - ! integer :: np !jpnote - integer(kind=timestepkind):: n integer :: i=0 @@ -756,18 +776,13 @@ subroutine integrate_gotm() #ifdef _ICE_ Qsw = I_0%value - !np = ubound(S,1) - call do_ice(h(nlev),dt,T(nlev),S(nlev),ta,precip%value,Qsw, & !h -- sum of dz ,, T -- temp of water !make sure arrays are not shaped for water + + call do_ice(h(nlev),dt,T(nlev),S(nlev),ta,precip%value,Qsw, & surface_fluxes,julianday,secondsofday,longitude, & latitude,I_0%value,airt%value,airp%value,hum%value, & - u10%value,v10%value,cloud%value,rho(nlev),rho_0,ql%method, & !sst,sss%value, !ql%method == longwave_radiation_method + u10%value,v10%value,cloud%value,rho(nlev),rho_0,ql%method, & hum_method,fluxes_method,albedo,heat%value) - !call do_ice(h(np),dt,T(np),S(np),ta,precip%value,Qsw, & !h -- sum of dz ,, T -- temp of water !make sure arrays are not shaped for water - ! surface_fluxes,julianday,secondsofday,longitude, & - ! latitude,I_0%value,airt%value,airp%value,hum%value, & - ! u10%value,v10%value,cloud%value,rho(np),rho_0,ql%method, & !sst,sss%value, !ql%method == longwave_radiation_method - ! hum_method,fluxes_method,albedo,heat%value) #endif ! reset some quantities @@ -805,7 +820,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) !jpnote passing ice_uvic_fs and ice_uvic_ff from ice + + 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 bb1d6e9a..62711e54 100644 --- a/src/gotm/register_all_variables.F90 +++ b/src/gotm/register_all_variables.F90 @@ -229,12 +229,12 @@ 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 +! 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 +! 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') @@ -258,14 +258,10 @@ subroutine register_stim_variables(nlev) 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') !changed from 1:nilay to 2:nilay+1 - - call fm%register('ice_uvic_zi', 'm', 'interface depth', standard_name='', dimensions=(/id_dim_zice/), data1d=ice_uvic_zi(2:nilay+2), category='ice') ! zi in old code corresponds to ice_uvic_zi in new code + 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') diff --git a/src/meanflow/salinity.F90 b/src/meanflow/salinity.F90 index c1db14f1..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: @@ -116,9 +116,9 @@ subroutine salinity(nlev,dt,cnpar,nus,gams,Fs,Ff) !jpnote passing in Fs and Ff f ! set boundary conditions DiffBcup = Neumann DiffBcdw = Neumann - !DiffSup = -S(nlev)*(precip%value+evap) !jpnote commenting + DiffSup = -S(nlev)*(precip%value+evap+Ff)-Fs !NSnote, check signs - DiffSup = -S(nlev)*(precip%value+evap+Ff)-Fs ! jpnote adding equation used in old code !need to pass in ff and fs for this + !DiffSup = -S(nlev)*(precip%value+evap) DiffSdw = _ZERO_