Merge PR #123 'ocean/time-varying-wind' into ocean/develop

This pull request adds time varying wind and pressure forcing using
Adrian Turner's framework implementation.

docs/ocean/design_docs/time-varying-wind.rst

@@ -0,0 +1,100 @@
+Setting up time-varying forcing for MPAS-O
+==========================================
+Implementation
+---------------
+ https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F 
+
+API
+---
+ https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/framework/mpas_forcing.F 
+
+Example datasets
+----------------
+ https://zenodo.org/record/1219191#.W9jWCKlRcUF with example from this dataset::
+
+    > pwd 
+    /Users/pwolfram/Downloads/MPAS-Seaice_test_dataset_V1/domain_QU240km
+
+    > ncdump -h LYq_six_hourly.2000.nc
+    netcdf LYq_six_hourly.2000 {
+    dimensions:
+        nCells = 1794 ;
+        StrLen = 64 ;
+         Time = UNLIMITED ; // (1460 currently)
+    variables:
+        char xtime(Time, StrLen) ;
+        double airTemperature(Time, nCells) ;
+        double airSpecificHumidity(Time, nCells) ;
+        double uAirVelocity(Time, nCells) ;
+        double vAirVelocity(Time, nCells) ;
+    }
+
+    > ncdump -v xtime LYq_six_hourly.2000.nc | tail
+      "2000-12-29_24:00:00                                             ",
+      "2000-12-30_06:00:00                                             ",
+      "2000-12-30_12:00:00                                             ",
+      "2000-12-30_18:00:00                                             ",
+      "2000-12-30_24:00:00                                             ",
+      "2000-12-31_06:00:00                                             ",
+      "2000-12-31_12:00:00                                             ",
+      "2000-12-31_18:00:00                                             ",
+      "2000-12-31_24:00:00                                             " ;
+    }
+
+Data structure
+--------------
+1. Global variables needed: `type (MPAS_forcing_group_type), pointer :: seaiceForcingGroups`
+
+2. Register against the `ForcingGroups` using the following initialization steps.  `ForcingGroups` is a pointer to a list of forcings, e.g., see  initialized group in https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L167 and members of the group https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L180 and https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L193 
+
+3. Streams file defines these source inputs, use example list of steps to get the data into a registry defined variable field / array.
+
+Initialization
+--------------
+Example of init steps (https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L232) 
+
+
+1. Create forcing group via `MPAS_forcing_init_group` https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L167 
+
+2. Initialize fields via `MPAS_forcing_init_field` https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L180 
+
+3. Initialize the data via reading the file with `MPAS_forcing_init_field_data`: https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L232 
+
+In time-step usage
+------------------
+Example of usage steps:
+
+1. `MPAS_forcing_get_forcing` Loop over all the individual forcing fields in the forcing group and get the data and perform the time interpolation. https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L389 
+
+2. `MPAS_forcing_get_forcing_time` Return the current forcing clock time for a forcing group. https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L517 
+
+3. `MPAS_forcing_write_restart_times` Loop over the forcing groups in the forcing group object and write out the forcing clock times to registry variables that are included in the restart stream. 'forcingGroupHead' is the forcing group object https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/shared/mpas_seaice_forcing.F#L2080 
+
+Registry changes needed (for restarts to work)
+----------------------------------------------
+Registry variables needed for restart:
+
+1. `nForcingGroupMax` https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/Registry.xml#L287 
+
+2. `nForcingGroupCounter` https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/Registry.xml#L3295 
+
+3. `forcingGroupNames` https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/Registry.xml#L3297  
+
+4. `forcingGroupRestartTimes` https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/Registry.xml#L3298 
+
+Register to restart streams:
+
+1. https://github.com/MPAS-Dev/MPAS-Model/blob/seaice/develop/src/core_seaice/Registry.xml#L1798
+
+ | `<var name="forcingGroupNames"/>` 
+ | `<var name="forcingGroupRestartTimes"/>`
+
+Additional resources
+--------------------
+M Maltrud PRs using this information: https://github.com/MPAS-Dev/MPAS/pulls?utf8=%E2%9C%93&q=is%3Apr+is%3Aclosed+author%3Amaltrud
+Old (depricated-- recommended to NOT USE) design doc: https://docs.google.com/document/d/1QtjmVCxLKS-S9Z_X8-WiHqXGn_yDu_VFHLA93ZZHg88/edit
+
+Time estimate (MPAS / MPAS-O novice)
+From start (this document) to merged PR:
+2 solid weeks
+4 weeks half-time

src/core_ocean/Registry.xml

@@ -713,6 +713,40 @@
 					possible_values="Any positive real number."
 		/>
 	</nml_record>
+	<nml_record name="time_varying_forcing" mode="init;forward">
+		<nml_option name="config_use_time_varying_atmospheric_forcing" type="logical" default_value=".false." units="unitless"
+			description="If true calculate input forcing fields."
+			possible_values=".true. or .false."
+		/>
+		<nml_option name="config_time_varying_atmospheric_forcing_type" type="character" default_value="WINDPRES" units="unitless"
+			description="Atmospheric forcing type."
+			possible_values="'WINDPRES'"
+		/>
+		<nml_option name="config_time_varying_atmospheric_forcing_start_time" type="character" default_value="0001-01-01_00:00:00" units="unitless"
+			description="Forcing time to use at the simulation start time"
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_atmospheric_forcing_reference_time" type="character" default_value="0001-01-01_00:00:00" units="unitless"
+			description="Reference time for the forcing"
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_atmospheric_forcing_cycle_start" type="character" default_value="0001-01-01_00:00:00" units="unitless"
+			description="Start time for the forcing cycle."
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_atmospheric_forcing_cycle_duration" type="character" default_value="2-00-00_00:00:00" units="unitless"
+			description="Duration of the forcing cycle."
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_atmospheric_forcing_interval" type="character" default_value="01:00:00" units="unitless"
+			description="Time between forcing inputs"
+			possible_values="'YYYY-MM-DD_HH:MM:SS"
+		/>
+		<nml_option name="config_time_varying_atmospheric_forcing_ramp" type="real" default_value="10.0" units="days"
+			description="Number of days to ramp up time varying forcing"
+			possible_values="Any positive real number"
+		/>
+	</nml_record>
 	<nml_record name="coupling" mode="init;forward">
 		<nml_option name="config_ssh_grad_relax_timescale" type="real" default_value="0.0" units="seconds"
 					description="Timescale for relaxation of the ssh gradient for coupling. A value of 0.0 (default) removes any relaxation and gives instantaneous response."
@@ -1207,6 +1241,7 @@
 		/>
 	</nml_record>
 	<packages>
+		<package name="timeVaryingAtmosphericForcingPKG" description="This package includes variables required for time varying atmospheric forcing"/>
 		<package name="variableShortwave" description="This package includes variables required to compute spatially variable shortwave extinction coefficients"/>
 
 		<package name="splitTimeIntegrator" description="This package includes variables required for either the split or unsplit explicit time integrators."/>
@@ -1602,6 +1637,17 @@
 			<var name="evaporationFlux"/>
 			<var name="rainFlux"/>
 		</stream>
+		<stream name="atmospheric_forcing"
+				type="input"
+				filename_template="atmospheric_forcing.nc"
+				filename_interval="none"
+				input_interval="none"
+				packages="timeVaryingAtmosphericForcingPKG"
+				immutable="true">
+			<var name="atmosPressure"/>
+			<var name="windSpeedU"/>
+			<var name="windSpeedV"/>
+		</stream>
 
 		<!-- Diagnostics streams for forward mode -->
 		<stream name="additional_output"
@@ -2977,6 +3023,24 @@
 			 packages="frazilIce"
 		/>
 	</var_struct>
+	<var_struct name="timeVaryingForcing" time_levs="1">
+		<var name="windSpeedU" type="real" dimensions="nCells Time" units="m/s"
+			description="Zonal (eastward) component of wind speed at cell centers from coupler. Positive easthward."
+			packages="timeVaryingAtmosphericForcingPKG"
+		/>
+		<var name="windSpeedV" type="real" dimensions="nCells Time" units="m/s"
+			description="Meridional (northward) component of wind speed at cell centers from coupler. Positive northward."
+			packages="timeVaryingAtmosphericForcingPKG"
+		/>
+		<var name="windSpeedMagnitude" type="real" dimensions="nCells Time" units="m/s"
+			description="Magnitude of wind speed at cell centers from coupler."
+			packages="timeVaryingAtmosphericForcingPKG"
+		/>
+		<var name="atmosPressure" type="real" dimensions="nCells Time" units="Pa"
+			description="Pressure at the sea surface due to the atmosphere."
+			packages="timeVaryingAtmosphericForcingPKG"
+		/>
+	</var_struct>
 	<var_struct name="scratch" time_levs="1">
 		<var name="normalThicknessFlux" persistence="scratch" type="real" dimensions="nVertLevels nEdges Time"
 			 units="m^{2} s^{-1}"

src/core_ocean/driver/mpas_ocn_core_interface.F

@@ -122,6 +122,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: inSituEOSActive
       logical, pointer :: variableShortwaveActive
       logical, pointer :: gmActive
+      logical, pointer :: timeVaryingAtmosphericForcingPKGActive
 
       type (mpas_pool_iterator_type) :: pkgItr
       logical, pointer :: packageActive
@@ -145,6 +146,8 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: config_use_freq_filtered_thickness
       logical, pointer :: config_use_frazil_ice_formation
       logical, pointer :: config_use_standardGM
+      logical, pointer :: config_use_time_varying_atmospheric_forcing
+
       character (len=StrKIND), pointer :: config_time_integrator
       character (len=StrKIND), pointer :: config_ocean_run_mode
       character (len=StrKIND), pointer :: config_pressure_gradient_type
@@ -288,6 +291,16 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
          gmActive = .true.
       end if
 
+      ! test for time-varying forcing
+      !
+      call mpas_pool_get_package(packagePool, 'timeVaryingAtmosphericForcingPKGActive', timeVaryingAtmosphericForcingPKGActive)
+      call mpas_pool_get_config(configPool, &
+                               'config_use_time_varying_atmospheric_forcing', &
+                                config_use_time_varying_atmospheric_forcing)
+      if (config_use_time_varying_atmospheric_forcing) then
+         timeVaryingAtmosphericForcingPKGActive = .true.
+      endif
+
       !
       ! call into analysis member driver to set analysis member packages
       !

src/core_ocean/mode_forward/mpas_ocn_forward_mode.F

@@ -76,6 +76,7 @@ module ocn_forward_mode
 
    use ocn_forcing
    use ocn_sea_ice
+   use ocn_time_varying_forcing
 
    use ocn_constants
 
@@ -484,6 +485,7 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
       logical, pointer :: config_use_ecosysTracers
       logical, pointer :: config_use_activeTracers_surface_restoring
       logical, pointer :: config_use_surface_salinity_monthly_restoring
+      logical, pointer :: config_use_time_varying_atmospheric_forcing
 
       ierr = 0
 
@@ -495,6 +497,8 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
                                                    config_use_activeTracers_surface_restoring)
       call mpas_pool_get_config(domain % configs, 'config_use_surface_salinity_monthly_restoring',  &
                                                    config_use_surface_salinity_monthly_restoring)
+      call MPAS_pool_get_config(domain % configs, 'config_use_time_varying_atmospheric_forcing', &
+                                                   config_use_time_varying_atmospheric_forcing)
 
       ! Eventually, dt should be domain specific
       timeStep = mpas_get_clock_timestep(domain % clock, ierr=ierr)
@@ -534,6 +538,10 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
         call mpas_timer_stop('io_monthly_surface_salinity')
       endif
 
+      ! initialize time-varying forcing
+      call ocn_time_varying_forcing_init(domain)
+      call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
+
       ! During integration, time level 1 stores the model state at the beginning of the
       !   time step, and time level 2 stores the state advanced dt in time by timestep(...)
       itimestep = 0
@@ -627,13 +635,14 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
                close(22)
             end if
 
-            if(trim(config_sw_absorption_type)=='ohlmann00') call ocn_shortwave_forcing_write_restart(domain)
-
-            if (config_use_ecosysTracers) call ocn_ecosys_forcing_write_restart(domain)
-
-            if (config_use_activeTracers_surface_restoring .and. config_use_surface_salinity_monthly_restoring)  &
-               call ocn_salinity_restoring_forcing_write_restart(domain)
+            ! write forcing restart if any type(s) time varying forcing is used
+            if (trim(config_sw_absorption_type)=='ohlmann00' .or. &
+                config_use_ecosysTracers .or. &
+                (config_use_activeTracers_surface_restoring .and. config_use_surface_salinity_monthly_restoring) .or. &
+                config_use_time_varying_atmospheric_forcing) then
 
+                call ocn_time_varying_forcing_write_restart_times(domain)
+            end if
          end if
 
          call mpas_timer_start('io_write')
@@ -674,6 +683,9 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
            endif
          endif
 
+         ! read in next time level data required for time-varying forcing
+         call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
+
          ! Validate that the state is OK to run with for the next timestep.
          call ocn_validate_state(domain, timeLevel=1)
 

src/core_ocean/shared/Makefile

@@ -64,7 +64,8 @@ OBJS = mpas_ocn_init_routines.o \
 	   mpas_ocn_forcing_restoring.o \
 	   mpas_ocn_time_average_coupled.o \
 	   mpas_ocn_sea_ice.o \
-	   mpas_ocn_framework_forcing.o
+	   mpas_ocn_framework_forcing.o \
+           mpas_ocn_time_varying_forcing.o
 
 all: $(OBJS)
 
@@ -194,6 +195,8 @@ mpas_ocn_time_average_coupled.o: mpas_ocn_constants.o
 
 mpas_ocn_framework_forcing.o:
 
+mpas_ocn_time_varying_forcing.o:  
+
 clean:
 	$(RM) *.o *.i *.mod *.f90
 

src/core_ocean/shared/mpas_ocn_constants.F

@@ -56,7 +56,8 @@ module ocn_constants
       latent_heat_fusion    ,&! lat heat of fusion (erg/g)
       latent_heat_fusion_mks,&! lat heat of fusion (J/kg)
       sea_ice_salinity      ,&! salinity of sea ice formed (psu)
-      ocn_ref_salinity        ! ocean reference salinity (psu)
+      ocn_ref_salinity      ,&! ocean reference salinity (psu)
+      atm_ref_pressure        ! standard sea level pressure (Pa)
 
    !  conversion factors
 
@@ -128,6 +129,7 @@ subroutine ocn_constants_init(configPool, packagePool)!{{{
        latent_heat_fusion_mks = 3.337e5_RKIND ! lat heat of fusion (J/kg)
        sea_ice_salinity   =  4.0_RKIND        ! (psu)
        ocn_ref_salinity   = 34.7_RKIND        ! (psu)
+       atm_ref_pressure   = 101325.0_RKIND    ! (Pa)
 
 
        !-----------------------------------------------------------------------