Merge pull request #189 from gustavo-marques/ncar_cfcs_implementation

* Implements a module to simulate CFCs via NUOPC cap

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -31,7 +31,7 @@ module MOM_cap_mod
 use MOM_ocean_model_nuopc,    only: ice_ocean_boundary_type
 use MOM_grid,                 only: ocean_grid_type, get_global_grid_size
 use MOM_ocean_model_nuopc,    only: ocean_model_restart, ocean_public_type, ocean_state_type
-use MOM_ocean_model_nuopc,    only: ocean_model_init_sfc
+use MOM_ocean_model_nuopc,    only: ocean_model_init_sfc, ocean_model_flux_init
 use MOM_ocean_model_nuopc,    only: ocean_model_init, update_ocean_model, ocean_model_end
 use MOM_ocean_model_nuopc,    only: get_ocean_grid, get_eps_omesh
 use MOM_cap_time,             only: AlarmInit
@@ -649,6 +649,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   ocean_public%is_ocean_pe = .true.
   call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(restartfiles))
 
+  ! GMM, this call is not needed for NCAR. Check with EMC.
+  ! If this can be deleted, perhaps we should also delete ocean_model_flux_init
+  call ocean_model_flux_init(ocean_state)
+
   call ocean_model_init_sfc(ocean_state, ocean_public)
 
   call mpp_get_compute_domain(ocean_public%domain, isc, iec, jsc, jec)
@@ -668,6 +672,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
              Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),&
              Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec),     &
              Ice_ocean_boundary% mi (isc:iec,jsc:jec),              &
+             Ice_ocean_boundary% ice_fraction (isc:iec,jsc:jec),    &
+             Ice_ocean_boundary% u10_sqr (isc:iec,jsc:jec),         &
              Ice_ocean_boundary% p (isc:iec,jsc:jec),               &
              Ice_ocean_boundary% lrunoff_hflx (isc:iec,jsc:jec),    &
              Ice_ocean_boundary% frunoff_hflx (isc:iec,jsc:jec),    &
@@ -689,6 +695,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   Ice_ocean_boundary%seaice_melt     = 0.0
   Ice_ocean_boundary%seaice_melt_heat= 0.0
   Ice_ocean_boundary%mi              = 0.0
+  Ice_ocean_boundary%ice_fraction    = 0.0
+  Ice_ocean_boundary%u10_sqr         = 0.0
   Ice_ocean_boundary%p               = 0.0
   Ice_ocean_boundary%lrunoff_hflx    = 0.0
   Ice_ocean_boundary%frunoff_hflx    = 0.0
@@ -747,6 +755,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "inst_pres_height_surface"   , "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl"                  , "will provide") !-> liquid runoff
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi"                  , "will provide") !-> ice runoff
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac"                   , "will provide") !-> ice fraction
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n"                  , "will provide") !-> wind^2 at 10m
   call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn"        , "will provide")
   !These are not currently used and changing requires a nuopc dictionary change

config_src/drivers/nuopc_cap/mom_cap_methods.F90

@@ -250,12 +250,27 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
   !----
   ! Note - preset values to 0, if field does not exist in importState, then will simply return
   ! and preset value will be used
-
   ice_ocean_boundary%mi(:,:) = 0._ESMF_KIND_R8
   call state_getimport(importState, 'mass_of_overlying_ice',  &
        isc, iec, jsc, jec, ice_ocean_boundary%mi,rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
+  !----
+  ! sea-ice fraction
+  !----
+  ice_ocean_boundary%ice_fraction(:,:) = 0._ESMF_KIND_R8
+  call state_getimport(importState, 'Si_ifrac',  &
+       isc, iec, jsc, jec, ice_ocean_boundary%ice_fraction, rc=rc)
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+  !----
+  ! 10m wind squared
+  !----
+  ice_ocean_boundary%u10_sqr(:,:) = 0._ESMF_KIND_R8
+  call state_getimport(importState, 'So_duu10n',  &
+       isc, iec, jsc, jec, ice_ocean_boundary%u10_sqr, rc=rc)
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
   !----
   ! Partitioned Stokes Drift Components
   !----

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -40,7 +40,6 @@ module MOM_ocean_model_nuopc
 use MOM_time_manager,        only : operator(/=), operator(<=), operator(>=)
 use MOM_time_manager,        only : operator(<), real_to_time_type, time_type_to_real
 use time_interp_external_mod,only : time_interp_external_init
-use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init
 use MOM_tracer_flow_control, only : call_tracer_flux_init
 use MOM_unit_scaling,        only : unit_scale_type
 use MOM_variables,           only : surface
@@ -240,14 +239,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
                                               !! tracer fluxes, and can be used to spawn related
                                               !! internal variables in the ice model.
   character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read
+
   ! Local variables
   real :: Rho0        ! The Boussinesq ocean density, in kg m-3.
   real :: G_Earth     ! The gravitational acceleration in m s-2.
   real :: HFrz        !< If HFrz > 0 (m), melt potential will be computed.
                       !! The actual depth over which melt potential is computed will
                       !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
                       !! If HFrz <= 0 (default), melt potential will not be computed.
-  logical :: use_melt_pot!< If true, allocate melt_potential array
+  logical :: use_melt_pot !< If true, allocate melt_potential array
+  logical :: use_CFC  !< If true, allocated arrays for surface CFCs.
+
 
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
@@ -366,10 +368,14 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
     use_melt_pot=.false.
   endif
 
+  call get_param(param_file, mdl, "USE_CFC_CAP", use_CFC, &
+                 default=.false., do_not_log=.true.)
+
   !   Consider using a run-time flag to determine whether to do the diagnostic
   ! vertical integrals, since the related 3-d sums are not negligible in cost.
   call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, &
-                              do_integrals=.true., gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot)
+                              do_integrals=.true., gas_fields_ocn=gas_fields_ocn, &
+                              use_meltpot=use_melt_pot, use_cfcs=use_CFC)
 
   call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, &
                             OS%forcing_CSp, OS%restore_salinity, OS%restore_temp)
@@ -416,6 +422,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
   endif
 
+  call extract_surface_state(OS%MOM_CSp, OS%sfc_state)
+
   call close_param_file(param_file)
   call diag_mediator_close_registration(OS%diag)
 

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -21,6 +21,7 @@ module MOM_surface_forcing_nuopc
 use MOM_forcing_type,     only : allocate_mech_forcing, deallocate_mech_forcing
 use MOM_get_input,        only : Get_MOM_Input, directories
 use MOM_grid,             only : ocean_grid_type
+use MOM_CFC_cap,          only : CFC_cap_fluxes
 use MOM_io,               only : slasher, write_version_number, MOM_read_data
 use MOM_restart,          only : register_restart_field, restart_init, MOM_restart_CS
 use MOM_restart,          only : restart_init_end, save_restart, restore_state
@@ -79,7 +80,7 @@ module MOM_surface_forcing_nuopc
                                 !! the correction for the atmospheric (and sea-ice)
                                 !! pressure limited by max_p_surf instead of the
                                 !! full atmospheric pressure.  The default is true.
-
+  logical :: use_CFC            !< enables the MOM_CFC_cap tracer package.
   real :: gust_const            !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa]
   logical :: read_gust_2d       !< If true, use a 2-dimensional gustiness supplied
                                 !! from an input file.
@@ -128,6 +129,7 @@ module MOM_surface_forcing_nuopc
 
   type(diag_ctrl), pointer :: diag                  !< structure to regulate diagnostic output timing
   character(len=200)       :: inputdir              !< directory where NetCDF input files are
+  character(len=200)       :: CFC_BC_file           !< filename with cfc11 and cfc12 data
   character(len=200)       :: salt_restore_file     !< filename for salt restoring data
   character(len=30)        :: salt_restore_var_name !< name of surface salinity in salt_restore_file
   logical                  :: mask_srestore         !< if true, apply a 2-dimensional mask to the surface
@@ -141,9 +143,13 @@ module MOM_surface_forcing_nuopc
                                                     !! temperature restoring fluxes. The masking file should be
                                                     !! in inputdir/temp_restore_mask.nc and the field should
                                                     !! be named 'mask'
+  character(len=30)        :: cfc11_var_name        !< name of cfc11 in CFC_BC_file
+  character(len=30)        :: cfc12_var_name        !< name of cfc11 in CFC_BC_file
   real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring
-  integer :: id_srestore = -1     !< id number for time_interp_external.
-  integer :: id_trestore = -1     !< id number for time_interp_external.
+  integer :: id_srestore = -1    !< id number for time_interp_external.
+  integer :: id_trestore = -1    !< id number for time_interp_external.
+  integer :: id_cfc11_atm = -1   !< id number for time_interp_external.
+  integer :: id_cfc12_atm = -1   !< id number for time_interp_external.
 
   ! Diagnostics handles
   type(forcing_diags), public :: handles
@@ -178,6 +184,8 @@ module MOM_surface_forcing_nuopc
   real, pointer, dimension(:,:) :: frunoff_hflx      =>NULL() !< heat content of frozen runoff [W/m2]
   real, pointer, dimension(:,:) :: p                 =>NULL() !< pressure of overlying ice and atmosphere
                                                               !< on ocean surface [Pa]
+  real, pointer, dimension(:,:) :: ice_fraction      =>NULL() !< mass of ice [nondim]
+  real, pointer, dimension(:,:) :: u10_sqr           =>NULL() !< wind speed squared at 10m [m2/s2]
   real, pointer, dimension(:,:) :: mi                =>NULL() !< mass of ice [kg/m2]
   real, pointer, dimension(:,:) :: ice_rigidity      =>NULL() !< rigidity of the sea ice, sea-ice and
                                                               !! ice-shelves, expressed as a coefficient
@@ -234,6 +242,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
   ! local variables
   real, dimension(SZI_(G),SZJ_(G)) :: &
+    cfc11_atm,     & !< CFC11 concentration in the atmopshere [???????]
+    cfc12_atm,     & !< CFC11 concentration in the atmopshere [???????]
     data_restore,  & !< The surface value toward which to restore [g/kg or degC]
     SST_anom,      & !< Instantaneous sea surface temperature anomalies from a target value [deg C]
     SSS_anom,      & !< Instantaneous sea surface salinity anomalies from a target value [g/kg]
@@ -293,7 +303,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   ! flux type has been used.
   if (fluxes%dt_buoy_accum < 0) then
     call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., &
-                               press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug)
+                               press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, &
+                               cfc=CS%use_CFC)
 
     call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed)
     call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed)
@@ -433,6 +444,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     enddo ; enddo
   endif
 
+
   ! Check that liquid runoff has a place to go
   if (CS%liquid_runoff_from_data .and. .not. associated(IOB%lrunoff)) then
     call MOM_error(FATAL, "liquid runoff is being added via data_override but "// &
@@ -496,6 +508,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
     if (associated(IOB%seaice_melt)) &
          fluxes%seaice_melt(i,j) = kg_m2_s_conversion * G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0)
 
+    ! sea ice fraction [nondim]
+    if (associated(IOB%ice_fraction)) &
+         fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0)
+
+    ! 10-m wind speed squared [m2/s2]
+    if (associated(IOB%u10_sqr)) &
+         fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0)
+
     fluxes%latent(i,j) = 0.0
     ! notice minus sign since fprec is positive into the ocean
     if (associated(IOB%fprec)) then
@@ -550,6 +570,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
      fluxes%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure.
   endif
 
+  ! CFCs
+  if (CS%use_CFC) then
+    call CFC_cap_fluxes(fluxes, sfc_state, G, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm)
+  endif
+
   if (associated(IOB%salt_flux)) then
     do j=js,je ; do i=is,ie
       fluxes%salt_flux(i,j)    = G%mask2dT(i,j)*(fluxes%salt_flux(i,j) + kg_m2_s_conversion*IOB%salt_flux(i-i0,j-j0))
@@ -1161,6 +1186,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
                  "coupler. This is used for testing and should be =1.0 for any "//&
                  "production runs.", default=1.0)
 
+  call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC, &
+                 default=.false., do_not_log=.true.)
+
   if (restore_salt) then
     call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, &
                  "The constant that relates the restoring surface fluxes to the relative "//&
@@ -1173,7 +1201,6 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
                  "The name of the surface salinity variable to read from "//&
                  "SALT_RESTORE_FILE for restoring salinity.", &
                  default="salt")
-
     call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", CS%salt_restore_as_sflux, &
                  "If true, the restoring of salinity is applied as a salt "//&
                  "flux instead of as a freshwater flux.", default=.false.)
@@ -1269,7 +1296,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
     enddo ; enddo
   endif
 
-  call time_interp_external_init
+  ! initialize time interpolator module
+  call time_interp_external_init()
 
 ! Optionally read a x-y gustiness field in place of a global
 ! constant.
@@ -1320,8 +1348,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
   call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, &
                  "If true, makes available diagnostics of fluxes from icebergs "//&
                  "as seen by MOM6.", default=.false.)
+
   call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles, &
-                                   use_berg_fluxes=iceberg_flux_diags)
+                                   use_berg_fluxes=iceberg_flux_diags, use_cfcs=CS%use_CFC)
 
   call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, &
                  "If true, allows flux adjustments to specified via the "//&
@@ -1355,6 +1384,27 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
     endif
   endif ; endif
 
+  ! Do not log these params here since they are logged in the CFC cap module
+  if (CS%use_CFC) then
+    call get_param(param_file, mdl, "CFC_BC_FILE", CS%CFC_BC_file, &
+                   "The file in which the CFC-11 and CFC-12 atm concentrations can be "//&
+                   "found (units must be parts per trillion), or an empty string for "//&
+                   "internal BC generation (TODO).", default=" ", do_not_log=.true.)
+    if ((len_trim(CS%CFC_BC_file) > 0) .and. (scan(CS%CFC_BC_file,'/') == 0)) then
+      ! Add the directory if CFC_BC_file is not already a complete path.
+      CS%CFC_BC_file = trim(slasher(CS%inputdir))//trim(CS%CFC_BC_file)
+      call get_param(param_file, mdl, "CFC11_VARIABLE", CS%cfc11_var_name, &
+                   "The name of the variable representing CFC-11 in  "//&
+                   "CFC_BC_FILE.", default="CFC_11", do_not_log=.true.)
+      call get_param(param_file, mdl, "CFC12_VARIABLE", CS%cfc12_var_name, &
+                   "The name of the variable representing CFC-12 in  "//&
+                   "CFC_BC_FILE.", default="CFC_12", do_not_log=.true.)
+
+      CS%id_cfc11_atm = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain)
+      CS%id_cfc12_atm = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain)
+    endif
+  endif
+
   ! Set up any restart fields associated with the forcing.
   call restart_init(param_file, CS%restart_CSp, "MOM_forcing.res")
   call restart_init_end(CS%restart_CSp)
@@ -1424,6 +1474,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   write(outunit,100) 'iobt%lrunoff        '   , mpp_chksum( iobt%lrunoff        )
   write(outunit,100) 'iobt%frunoff        '   , mpp_chksum( iobt%frunoff        )
   write(outunit,100) 'iobt%p              '   , mpp_chksum( iobt%p              )
+  write(outunit,100) 'iobt%ice_fraction   '   , mpp_chksum( iobt%ice_fraction   )
+  write(outunit,100) 'iobt%u10_sqr        '   , mpp_chksum( iobt%u10_sqr        )
   if (associated(iobt%ustar_berg)) &
     write(outunit,100) 'iobt%ustar_berg     ' , mpp_chksum( iobt%ustar_berg )
   if (associated(iobt%area_berg)) &