Uvic fabm jp

.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:

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)

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

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')       

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_