Uvic ice dev

.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

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,18 +1018,17 @@ 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)
 #else
       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))

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

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,15 +61,21 @@ subroutine longwave_radiation(method,dlat,tw,ta,cloud,ql)
 
    REALTYPE                  :: ccf
    REALTYPE                  :: x1,x2,x3
-!
 !EOP
 !-----------------------------------------------------------------------
 !BOC
 !  calculate cloud correction factor,fortran counts from 1 !
    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

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