! $Id: GEOS_TurbulenceGridComp.F90,v 1.88.2.3.2.1.2.2.2.6 2007/10/03 19:44:02 f4mjs Exp $
#include "MAPL_Generic.h"
!=============================================================================
module GEOS_TurbulenceGridCompMod 1,1
!BOP
! !MODULE: GEOS_Turbulence --- An GEOS generic atmospheric turbulence component
! !USES:
use ESMF_Mod
use GEOS_Mod
use LockEntrain
implicit none
private
! !PUBLIC MEMBER FUNCTIONS:
public SetServices
! !DESCRIPTION:
!
! {\tt GEOS\_TurbulenceGridComp} computes atmospheric tendencies due to turbulence.
! Its physics is a combination of the first-order scheme of Louis---for stable PBLs
! and free atmospheric turbulence---with a modified version of the non-local-K
! scheme proposed by Lock for unstable and cloud-topped boundary layers.
! In addition to diffusive tendencies, it adds the effects orographic form drag
! for features with horizontal scales of 2 to 20 km following Beljaars et al. (2003,
! ECMWF Tech. Memo. 427).
!
!\vspace{12 pt}
!\noindent
!{\bf Grid Considerations}
!
! Like all GEOS\_Generic-based components, it works on an inherited
! 3-dimensional ESMF grid. It assumes that the first two (inner) dimensions span the
! horizontal and the third (outer) dimension is the vertical. In the horizontal,
! one or both dimensions can be degenerate, effectively supporting
! single-columns (1-D), and slices (2-D). No horizontal dimension needs to be
! aligned with a particular coordinate. In the vertical, the only assumption
! is that columns are indexed from top to bottom.
!
!\vspace{12 pt}
!\noindent
!{\bf Methods}
!
! {\tt GEOS\_TurbulenceGridComp} uses the default Initialize and Finalize methods
! of GEOS\_Generic. It has a 2-stage Run method that can be used in conjunction with
! two-stage surface calculations to implement semi-implicit time differencing.
!
!\vspace{12 pt}
!\noindent
!{\bf Time Behavior}
!
! {\tt GEOS\_TurbulenceGridComp} assumes both run stages will be invoked every
! RUN\_DT seconds, where RUN\_DT is required in the configuration. On this interval
! both run stages will perform diffusion updates using diffusivities found in the
! internal state. The diffusivities in the internal state may be refreshed intermitently
! by specifying MY\_STEP and ACCUMINT in the configuration. Accumulated imports used
! in the intermittent refreshing are valid only on MY\_STEP intervals. Currently the
! origin of these intervals is the beginning of the run. Accumulation of these imports
! is done for a period ACCUMINT prior to the valid time. Both ACCUMINT and MY\_STEP are
! in seconds.
!
!\vspace{12 pt}
!\noindent
!{\bf Working with Bundles and Friendlies}
!
! {\tt GEOS\_TurbulenceGridComp} works on bundles of quantities to be diffused
! and with corresponding bundles of their tendencies, surface values, etc.
! These bundles may contain an arbitrary number of conservative quantities and
! no requirements or restrictions are placed on what quantities they contain.
! Quantities required for the calculation, such as pressures, stability, etc
! are passed separately from the diffused quantities. Little distinction is made
! of what is in the bundle, except that needed to decide what diffusivity applies
! to the quantity and in what form its effects are implemented.
!
! Quantities to be diffused can be marked as "Friendly-for-diffusion". In that case,
! {\tt GEOS\_TurbulenceGridComp} directly updates the quantity; otherwise it
! merely computes its tendency, placing it in the appropriate bundle and treating
! the quantity itself as read-only.
!
! In working with bundled quantities, corresponding fields must appear in the
! same order in all bundles. Some of these fields, however,
! may be ``empty'' in the sense that the data pointer has not been allocated.
!
! {\tt GEOS\_TurbulenceGridComp} works with six bundles; three in the import
! state and three in the export state. The import bundles are:
! \begin{itemize}
! \item[]
! \makebox[1in][l]{\bf TR}
! \parbox[t]{4in}{The quantity being diffused.}
! \item[]
! \makebox[1in][l]{\bf TRG}
! \parbox[t]{4in}{The surface (ground) value of the quantity being diffused.
! (Used only by Run2)}
! \item[]
! \makebox[1in][l]{\bf DTG}
! \parbox[t]{4in}{The change of TRG during the time step. (Used only by Run2)}
! \end{itemize}
!
! The export bundles are:
! \begin{itemize}
! \item[]
! \makebox[1in][l]{\bf TRI}
! \parbox[t]{4in}{The tendency of the quantity being diffused.
! (Produced by Run1, updated by Run2.) }
! \item[]
! \makebox[1in][l]{\bf FSTAR}
! \parbox[t]{4in}{After Run1, the ``preliminary'' (i.e., at the original surface
! value) surface flux of the diffused quantity; after Run2, its final value.
! (Produced by Run1, updated by Run2)}
! \item[]
! \makebox[1in][l]{\bf DFSTAR}
! \parbox[t]{4in}{The change of preliminary FSTAR per unit change in the
! surface value. (Produced by Run1)}
! \end{itemize}
!
! All fields in the export bundles are checked for associated pointers before being
! updated.
!
! Fields in the TR bundle can have four attributes:
! \begin{itemize}
! \item FriendlyTo[{\it Component Name}]: default=false --- If true, TR field is updated.
! \item WeightedTendency: default=true --- If true, tendencies (TRI) are pressure-weighted
! \item DiffuseLike: ('S','Q','M') default='S' --- Use mixing coefficients for either
! heat, moisture or momentum.
! \end{itemize}
!
! Only fields in the TR bundle are checked for friendly status. Non-friendly
! fields in TR and all other bundles are treated with the usual Import/Export
! rules.
!
!\vspace{12 pt}
!\noindent
!{\bf Other imports and exports}
!
! In addition to the updates of these bundles, {\tt GEOS\_TurbulenceGridComp} produces
! a number of diagnostic exports, as well as frictional heating contributions. The latter
! are NOT added by {\tt GEOS\_TurbulenceGridComp}, but merely exported to be added
! elsewhere in the GCM.
!
!\vspace{12 pt}
!\noindent
!{\bf Two-Stage Interactions with the Surface}
!
! The two-stage scheme for interacting with the surface module is as follows:
! \begin{itemize}
! \item The first run stage takes the surface values of the diffused quantities
! and the surface exchange coefficients as input. These are, of course, on the
! grid turbulence is working on.
! \item It then does the full diffusion calculation assuming the surface values are
! fixed, i.e., the explicit surface case. In addition, it also computes derivatives of the
! tendencies wrt surface values. These are to be used in the second stage.
! \item The second run stage takes the increments of the surface values as inputs
! and produces the final results, adding the implicit surface contributions.
! \item It also computes the frictional heating due to both implicit and explicit
! surface contributions.
! \end{itemize}
!
!\vspace{12 pt}
!\noindent
!{\bf GEOS-5 Specific Aspects}
!
! In GEOS-5, {\tt GEOS\_TurbulenceGridComp} works on the atmosphere's lat-lon grid,
! while surface quantities are computed during the first run stage of the each of
! the tiled surface components. The tiled quantities are properly aggregated to
! the GEOS-5 lat-lon grid by the first stage of {\tt GEOS\_SurfaceGridComp}, which
! is called immediately before the first run stage of {\tt GEOS\_TurbulenceGridComp}.
!
!EOP
contains
!=============================================================================
!=============================================================================
!=============================================================================
!=============================================================================
!BOP
! !IROUTINE: SetServices -- Sets ESMF services for this component
! !DESCRIPTION: This version uses the {\tt GEOS\_GenericSetServices}, which sets
! the Initialize and Finalize services to generic versions. It also
! allocates our instance of a generic state and puts it in the
! gridded component (GC). Here we only set the two-stage run method and
! declare the data services.
! \newline
! !REVISION HISTORY:
! ??Jul2006 E.Novak./Todling - Added output defining TLM/ADM trajectory
! !INTERFACE:
subroutine SetServices ( GC, RC ),1759
! !ARGUMENTS:
type(ESMF_GridComp), intent(INOUT) :: GC ! gridded component
integer, optional :: RC ! return code
!EOP
!=============================================================================
!
! ErrLog Variables
character(len=ESMF_MAXSTR) :: IAm
integer :: STATUS
character(len=ESMF_MAXSTR) :: COMP_NAME
! Local derived type aliases
type (ESMF_Config ) :: CF
! Locals
integer :: MY_STEP
integer :: ACCUMINT
real :: DT
!=============================================================================
! Begin...
! Get my name and set-up traceback handle
! ---------------------------------------
Iam = 'SetServices'
call ESMF_GridCompGet( GC, NAME=COMP_NAME, RC=STATUS )
VERIFY_(STATUS)
Iam = trim(COMP_NAME) // Iam
! Set the Run entry points
! ------------------------
call MAPL_GridCompSetEntryPoint ( GC, ESMF_SETRUN, Run1, RC=STATUS )
VERIFY_(STATUS)
call MAPL_GridCompSetEntryPoint ( GC, ESMF_SETRUN, Run2, RC=STATUS )
VERIFY_(STATUS)
! Get the intervals
! -----------------
call ESMF_GridCompGet( GC, CONFIG = CF, RC=STATUS )
VERIFY_(STATUS)
call ESMF_ConfigGetAttribute( CF, DT, Label="RUN_DT:" , RC=STATUS)
VERIFY_(STATUS)
call ESMF_ConfigGetAttribute( CF, DT, Label=trim(COMP_NAME)//"_DT:" , default=DT, RC=STATUS)
MY_STEP = nint(DT)
call ESMF_ConfigGetAttribute( CF, DT, Label=trim(COMP_NAME)//"Avrg:", default=DT, RC=STATUS)
ACCUMINT = nint(DT)
! Set the state variable specs.
! -----------------------------
!BOP
! !IMPORT STATE:
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'PLE', &
LONG_NAME = 'air_pressure', &
UNITS = 'Pa', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'T', &
LONG_NAME = 'air_temperature', &
UNITS = 'K', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'TH', &
LONG_NAME = 'potential_temperature', &
UNITS = 'K', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'QV', &
LONG_NAME = 'specific_humidity', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'QLLS', &
LONG_NAME = 'liquid_condensate_mixing_ratio', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'QILS', &
LONG_NAME = 'frozen_condensate_mixing_ratio', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'CLLS', &
LONG_NAME = 'cloud_fraction', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'QLCN', &
LONG_NAME = 'liquid_condensate_mixing_ratio', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'QICN', &
LONG_NAME = 'frozen_condensate_mixing_ratio', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'CLCN', &
LONG_NAME = 'cloud_fraction', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'U', &
LONG_NAME = 'eastward_wind', &
UNITS = 'm s-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'V', &
LONG_NAME = 'northward_wind', &
UNITS = 'm s-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'CT', &
LONG_NAME = 'surface_heat_exchange_coefficient', &
UNITS = 'kg m-2 s-1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'CQ', &
LONG_NAME = 'surface_moisture_exchange_coefficient', &
UNITS = 'kg m-2 s-1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'CM', &
LONG_NAME = 'surface_momentum_exchange_coefficient', &
UNITS = 'kg m-2 s-1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'BSTAR', &
LONG_NAME = 'surface_bouyancy_scale', &
UNITS = 'm s-2', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'USTAR', &
LONG_NAME = 'surface_velocity_scale', &
UNITS = 'm s-1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'FRLAND', &
LONG_NAME = 'land_fraction', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'RADLW', &
LONG_NAME = 'air_temperature_tendency_due_to_longwave',&
UNITS = 'K s-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'RADLWC', &
LONG_NAME = 'clearsky_air_temperature_tendency_lw',&
UNITS = 'K s-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
AVERAGING_INTERVAL = ACCUMINT, &
REFRESH_INTERVAL = MY_STEP, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'PREF', &
LONG_NAME = 'reference_air_pressure', &
UNITS = 'Pa', &
DIMS = MAPL_DimsVertOnly, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'VARFLT', &
LONG_NAME = 'variance_of_filtered_topography', &
UNITS = 'm+2', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'TR', &
LONG_NAME = 'diffused_quantities', &
UNITS = 'X', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_BundleItem, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'TRG', &
LONG_NAME = 'surface_values_of_diffused_quantity',&
UNITS = 'X', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
DATATYPE = MAPL_BundleItem, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddImportSpec(GC, &
SHORT_NAME = 'DTG', &
LONG_NAME = 'change_of_surface_values_of_diffused_quantity',&
UNITS = 'X', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
DATATYPE = MAPL_BundleItem, &
RC=STATUS )
VERIFY_(STATUS)
! !EXPORT STATE:
call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'TRI', &
LONG_NAME = 'diffusion_tendencies', &
UNITS = 'X kg m-2 s-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_BundleItem, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'FSTAR', &
LONG_NAME = 'surface_fluxes', &
UNITS = 'X kg m-2 s-1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
DATATYPE = MAPL_BundleItem, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'DFSTAR', &
LONG_NAME = 'change_of_surface_fluxes_for_unit_change_of_surface_value',&
UNITS = 'kg m-2 s-1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
DATATYPE = MAPL_BundleItem, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'air_temperature', &
UNITS = 'K', &
SHORT_NAME = 'T', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'eastward_wind', &
UNITS = 'm/s', &
SHORT_NAME = 'U', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'northward_wind', &
UNITS = 'm/s', &
SHORT_NAME = 'V', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'specific_humidity', &
UNITS = '1', &
SHORT_NAME = 'QV', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'total_momentum_diffusivity', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'KM', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'total_scalar_diffusivity', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'KH', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'Richardson_number_from_Louis', &
UNITS = '1', &
SHORT_NAME = 'RI', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'bulk_shear_from_Louis', &
UNITS = 's-1', &
SHORT_NAME = 'DU', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'scalar_diffusivity_from_Louis', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'KHLS', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'momentum_diffusivity_from_Louis', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'KMLS', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'surface_driven_scalar_diffusivity_from_Lock_scheme', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'KHSFC', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'radiation_driven_scalar_diffusivity_from_Lock_scheme', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'KHRAD', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'cloudy_LW_radiation_tendency_used_by_Lock_scheme', &
UNITS = 'K s-1', &
SHORT_NAME = 'LWCRT', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'entrainment_heat_diffusivity_from_Lock', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'EKH', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'entrainment_momentum_diffusivity_from_Lock', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'EKM', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'Blackadar_length_scale_for_scalars', &
UNITS = 'm', &
SHORT_NAME = 'ALH', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationEdge, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'p-weighted_frictional_heating_rate_from_diffusion', &
UNITS = 'K s-1 Pa', &
SHORT_NAME = 'INTDIS', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'p-weighted_frictional_heating_rate_from_orographic_drag',&
UNITS = 'K s-1 Pa', &
SHORT_NAME = 'TOPDIS', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'p-weighted_frictional_heating_rate_from_surface_drag', &
UNITS = 'K s-1 Pa', &
SHORT_NAME = 'SRFDIS', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec ( gc, &
SHORT_NAME = 'KETRB', &
LONG_NAME = 'vertically_integrated_kinetic_energy_tendency_across_turbulence',&
UNITS = 'W m-2', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec ( gc, &
SHORT_NAME = 'KESRF', &
LONG_NAME = 'vertically_integrated_kinetic_energy_dissipation_due_to_surface_friction',&
UNITS = 'W m-2', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec ( gc, &
SHORT_NAME = 'KEINT', &
LONG_NAME = 'vertically_integrated_kinetic_energy_dissipation_due_to_diffusion',&
UNITS = 'W m-2', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec ( gc, &
SHORT_NAME = 'KETOP', &
LONG_NAME = 'vertically_integrated_kinetic_energy_dissipation_due_to_topographic_friction',&
UNITS = 'W m-2', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'entrainment_velocity_from_surface_plume', &
UNITS = 'm s-1', &
SHORT_NAME = 'WESFC', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'entrainment_velocity_from_radiation', &
UNITS = 'm s-1', &
SHORT_NAME = 'WERAD', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'entrainment_velocity_from_buoy_rev', &
UNITS = 'm s-1', &
SHORT_NAME = 'WEBRV', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'Buoyancy_jump_across_inversion', &
UNITS = 'm s-2', &
SHORT_NAME = 'DBUOY', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'turbulent_velocity_scale_for_sfc', &
UNITS = 'm s-1', &
SHORT_NAME = 'VSCSFC', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'turbulent_velocity_scale_for_cooling', &
UNITS = 'm s-1', &
SHORT_NAME = 'VSCRAD', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'turbulent_velocity_scale_for_buoy_rev', &
UNITS = 'm s-1', &
SHORT_NAME = 'VSCBRV', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'turbulent_entrainment_diff_from_cooling', &
UNITS = 'm+2 s-1', &
SHORT_NAME = 'KERAD', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'cloud_top_radiative_forcing', &
UNITS = 'W m-2', &
SHORT_NAME = 'CLDRF', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'pbltop_height_as_diagnosed_from_KH', &
UNITS = 'm', &
SHORT_NAME = 'ZPBL', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'pbltop_pressure', &
UNITS = 'Pa', &
SHORT_NAME = 'PPBL', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'pbltop_height_for_sfc_plume_LOCK', &
UNITS = 'm', &
SHORT_NAME = 'ZSML', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'depth_for_rad/brv_plume_LOCK', &
UNITS = 'm', &
SHORT_NAME = 'ZRADML', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'hght_of_base_for_rad/brv_plume_LOCK', &
UNITS = 'm', &
SHORT_NAME = 'ZRADBS', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'pbltop_cloud_depth_LOCK', &
UNITS = 'm', &
SHORT_NAME = 'ZCLD', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'pbltop_cloud_top_height_LOCK', &
UNITS = 'm', &
SHORT_NAME = 'ZCLDTOP', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'optimal_mixture_fraction_for_BRV', &
UNITS = '1', &
SHORT_NAME = 'CHIS', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 's_of_optimal_mixture_for_BRV', &
UNITS = 'J kg-1', &
SHORT_NAME = 'SMIXT', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'Scaled_Del_s_at_Cloud_top', &
UNITS = 'K', &
SHORT_NAME = 'DELSINV', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'Siems_buoy_rev_parameter', &
UNITS = '1', &
SHORT_NAME = 'DSIEMS', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'Return_codes_for_Lock_top_driven_plume', &
UNITS = '1', &
SHORT_NAME = 'RADRCODE', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'matrix_diagonal_ak_for_scalars_over_dt', &
SHORT_NAME = 'AKSODT', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'matrix_diagonal_ck_for_scalars_over_dt', &
SHORT_NAME = 'CKSODT', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'matrix_diagonal_ak_for_moisture_over_dt', &
SHORT_NAME = 'AKQODT', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'matrix_diagonal_ck_for_moisture_over_dt', &
SHORT_NAME = 'CKQODT', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'matrix_diagonal_ak_for_winds_over_dt', &
SHORT_NAME = 'AKVODT', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddExportSpec(GC, &
LONG_NAME = 'matrix_diagonal_ck_for_winds_over_dt', &
SHORT_NAME = 'CKVODT', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
RC=STATUS )
VERIFY_(STATUS)
! !INTERNAL STATE:
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_ahat_for_scalars', &
SHORT_NAME = 'AKS', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_bhat_for_scalars', &
SHORT_NAME = 'BKS', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_c_for_scalars', &
SHORT_NAME = 'CKS', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'sensitivity_of_tendency_to_surface_value_for_scalars', &
SHORT_NAME = 'DKS', &
UNITS = 's-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_ahat_for_moisture', &
SHORT_NAME = 'AKQ', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_bhat_for_moisture', &
SHORT_NAME = 'BKQ', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_c_for_moisture', &
SHORT_NAME = 'CKQ', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'sensitivity_of_tendency_to_surface_value_for_moisture', &
SHORT_NAME = 'DKQ', &
UNITS = 's-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_ahat_for_winds', &
SHORT_NAME = 'AKV', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_bhat_for_winds', &
SHORT_NAME = 'BKV', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'matrix_diagonal_c_for_winds', &
SHORT_NAME = 'CKV', &
UNITS = '1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'sensitivity_of_tendency_to_surface_value_for_winds', &
SHORT_NAME = 'DKV', &
UNITS = 's-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'momentum_mixing_factor', &
SHORT_NAME = 'EKV', &
UNITS = 'Pa s-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
call MAPL_AddInternalSpec(GC, &
LONG_NAME = 'topographic_roughness_factor', &
SHORT_NAME = 'FKV', &
UNITS = 'Pa s-1', &
DIMS = MAPL_DimsHorzVert, &
VLOCATION = MAPL_VLocationCenter, &
DATATYPE = MAPL_NoRestart, &
RC=STATUS )
VERIFY_(STATUS)
!EOP
! Set the Profiling timers
! ------------------------
call MAPL_TimerAdd(GC, name="-RUN1" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="--DIFFUSE" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="--REFRESHKS" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="---PRELIMS" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="---LOUIS" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="---LOCK" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="---BELJAARS" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="---DECOMP" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="-RUN2" ,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerAdd(GC, name="--UPDATE" ,RC=STATUS)
VERIFY_(STATUS)
! Set generic init and final methods
! ----------------------------------
call MAPL_GenericSetServices ( GC, RC=STATUS)
VERIFY_(STATUS)
RETURN_(ESMF_SUCCESS)
end subroutine SetServices
!=============================================================================
!=============================================================================
!=============================================================================
!=============================================================================
!=============================================================================
!BOP
! !IROUTINE: RUN1 -- First run stage for the {\tt MAPL_TurbulenceGridComp} component
! !INTERFACE:
subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC ),83
! !ARGUMENTS:
type(ESMF_GridComp), intent(inout) :: GC
type(ESMF_State), intent(inout) :: IMPORT
type(ESMF_State), intent(inout) :: EXPORT
type(ESMF_Clock), intent(inout) :: CLOCK
integer, optional, intent( out) :: RC
! !DESCRIPTION: The first run stage of {\tt GEOS\_TurbulenceGridComp} computes the diffusivities,
! sets-up the matrix for a backward-implicit computation of the surface fluxes,
! and solves this system for a fixed surface value of the diffused quantity. Run1
! takes as inputs the surface exchange coefficients (i.e., $\rho |U| C_{m,h,q}$) for
! momentun, heat, and moisture, as well as the pressure, temperature, moisture,
! and winds for the sounding. These are used only for computing the diffusivities
! and, as explained above, are not the temperatures, moistures, etc. being diffused.
!
! The computation of turbulence fluxes for fixed surface values is done at every
! time step in the contained subroutine {\tt DIFFUSE}; but the computation of
! diffusivities and orographic drag coefficients, as well as the set-up of the
! vertical difference matrix and its LU decomposition
! can be done intermittently for economy in the contained subroutine {\tt REFRESH}.
! The results of this calculation are stored in an internal state.
! Run1 also computes the sensitivity of the
! atmospheric tendencies and the surface flux to changes in the surface value.
!
! The diffusivities are computed by calls to {\tt LOUIS\_KS} and {\tt ENTRAIN}, which
! compute the Louis et al. (1983) and Lock (2000) diffusivities. The Louis
! diffusivities are computed for all conditions, and {\tt ENTRAIN} overrides them
! where appropriate. Lock can be turned off from the resource file.
!
!EOP
! ErrLog Variables
character(len=ESMF_MAXSTR) :: IAm
integer :: STATUS
character(len=ESMF_MAXSTR) :: COMP_NAME
! Local derived type aliases
type (MAPL_MetaComp), pointer :: MAPL
type (ESMF_Config ) :: CF
type (ESMF_State ) :: INTERNAL
type (ESMF_Alarm ) :: ALARM
! Local variables
real, dimension(:,:,:), pointer :: AKS, BKS, CKS, DKS
real, dimension(:,:,:), pointer :: AKQ, BKQ, CKQ, DKQ
real, dimension(:,:,:), pointer :: AKV, BKV, CKV, DKV, EKV, FKV
real, dimension(:,:,:), pointer :: PLE
real, dimension(:,: ), pointer :: CU, CT, CQ
integer :: IM, JM, LM
real :: DT
! Begin...
!---------
! Get my name and set-up traceback handle
! ---------------------------------------
call ESMF_GridCompGet( GC, NAME=COMP_NAME, RC=STATUS )
VERIFY_(STATUS)
Iam = trim(COMP_NAME) // 'Run1'
! Get my internal MAPL_Generic state
!-----------------------------------
call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerOn(MAPL,"TOTAL")
call MAPL_TimerOn(MAPL,"-RUN1")
! Get parameters from generic state.
!-----------------------------------
call MAPL_Get(MAPL, &
IM=IM, JM=JM, LM=LM, &
RUNALARM=ALARM, &
INTERNAL_ESMF_STATE=INTERNAL, &
RC=STATUS )
VERIFY_(STATUS)
! Get configuration from component
!---------------------------------
call ESMF_GridCompGet( GC, CONFIG = CF, RC=STATUS )
VERIFY_(STATUS)
! Get all pointers that are needed by both REFRESH and DIFFUSE
!-------------------------------------------------------------
! Get pressure structure; this is instantaneous.
!-----------------------------------------------
call MAPL_GetPointer(IMPORT, PLE, 'PLE', RC=STATUS)
VERIFY_(STATUS)
! Get surface exchange coefficients
!----------------------------------
call MAPL_GetPointer(IMPORT, CU, 'CM', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, CT, 'CT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, CQ, 'CQ', RC=STATUS)
VERIFY_(STATUS)
! Get pointers from internal state
!---------------------------------
call MAPL_GetPointer(INTERNAL, AKS, 'AKS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, BKS, 'BKS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, CKS, 'CKS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, DKS, 'DKS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, AKQ, 'AKQ', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, BKQ, 'BKQ', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, CKQ, 'CKQ', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, DKQ, 'DKQ', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, AKV, 'AKV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, BKV, 'BKV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, CKV, 'CKV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, DKV, 'DKV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, EKV, 'EKV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, FKV, 'FKV', RC=STATUS)
VERIFY_(STATUS)
! Get application's timestep from configuration
!----------------------------------------------
call ESMF_ConfigGetAttribute(CF, DT, Label="RUN_DT:" , RC=STATUS)
VERIFY_(STATUS)
! If its time, do the refresh
! ---------------------------
if ( ESMF_AlarmIsRinging(ALARM, rc=status) ) then
VERIFY_(STATUS)
call ESMF_AlarmRingerOff(ALARM, RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerOn (MAPL,"--REFRESHKS")
call REFRESH(IM,JM,LM,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerOff(MAPL,"--REFRESHKS")
endif
! Solve the free atmosphere problem
! ---------------------------------
call MAPL_TimerOn (MAPL,"--DIFFUSE")
call DIFFUSE(IM,JM,LM,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerOff(MAPL,"--DIFFUSE")
! All done with RUN1
!--------------------
call MAPL_TimerOff(MAPL,"-RUN1")
call MAPL_TimerOff(MAPL,"TOTAL")
RETURN_(ESMF_SUCCESS)
contains
!=============================================================================
!=============================================================================
!BOP
! !CROUTINE: REFRESH -- Refreshes diffusivities.
! !INTERFACE:
subroutine REFRESH(IM,JM,LM,RC) 1,18
! !ARGUMENTS:
integer, intent(IN) :: IM,JM,LM
integer, optional, intent(OUT) :: RC
! !DESCRIPTION:
! {\tt REFRESH} can be called intermittently to compute new values of the
! diffusivities. In addition it does all possible calculations that depend
! only on these. In particular, it sets up the semi-implicit tridiagonal
! solver in the vertical and does the LU decomposition. It also includes the
! local effects of orographic drag, so that it to is done implicitly.
!
! Diffusivities are first computed with the Louis scheme ({\tt LOUIS\_KS}),
! and then, where appropriate,
! they are overridden by the Lock values ({\tt ENTRAIN}).
! Once diffusivities are computed, {\tt REFRESH} sets-up the tridiagonal
! matrices for the semi-implicit vertical diffusion calculation and performs
! their $LU$ decomposition.
!
! {\tt REFRESH} requires surface exchange coefficients for heat, moisture, and
! momentum, The calculations in the interior are also
! done for momentum, heat, and water diffusion. Heat and water mixing
! coefficients differ only at the surface, but these affect the entire $LU$
! decomposition, and so all three decompositions are saved in the internal state.
!
! For a conservatively diffused quantity $q$, we have
! $$
! \frac{\partial q}{\partial t} = -g \frac{\partial }{\partial p}
! \left(\rho K_q \frac{\partial q}{\partial z} \right)
! $$
! In finite difference form, using backward time differencing, this becomes
! $$
! \begin{array}{rcl}
! {q^{n+1}_l - q^{n}_l} & = & - \frac{g}{\delta_l p}^*
! \delta_l \left[
! \left( \frac{\Delta t \rho K_q}{\delta_l z} \right)^* (\delta_l q)^{n+1} \right] \\
! &&\\
! & = & - \alpha_l ( \beta_{l+\frac{1}{2}}(q_{l+1}-q_l)^{n+1} -
! \beta_{l-\frac{1}{2}}(q_l-q_{l-1})^{n+1} ) \\
! &&\\
! \alpha_l & = & \frac{g \Delta t}{(p_{l+\frac{1}{2}}-p_{l-\frac{1}{2}})^*} \\
! &&\\
! \beta_{l+\frac{1}{2}} & = & \left( \frac{ (\rho K_q)^*_{l+\frac{1}{2}}}{(z_{l+1}-z_{l})^*} \right) \\
! \end{array}
! $$
! where the subscripts denote levels, superscripts denote times, and the $*$ superscript
! denotes evaluation at the refresh time.
! The following tridiagonal set is then solved for $q^{n+1}_l$:
! $$
! a_l q_{l-1} + b_l q_l + c_l q_{l+1} = q_l
! $$
! where
! $$
! \begin{array}{rcl}
! a_l & = & \alpha_l \beta_{l-\frac{1}{2}} \\
! c_l & = & \alpha_l \beta_{l+\frac{1}{2}} \\
! b_l & = & 1 - a_l - c_l.
! \end{array}
! $$
! At the top boundary, we assume $K_q=0$, so $ \beta_{\frac{1}{2}}=0$ and $a_1=0$.
! At the surface, $ \beta_{L+\frac{1}{2}}= \rho_s |U|_s C_{m,h,q}$, the surface exchange coefficient.
!
!EOP
character(len=ESMF_MAXSTR) :: IAm='Refresh'
integer :: STATUS
real, dimension(:,:,:), pointer :: TH, U, V, Q, T, RI, DU, RADLW, RADLWC, LWCRT
real, dimension(:,: ), pointer :: VARFLT
real, dimension(:,:,:), pointer :: KH, KM, QLLS, QILS, CLLS, QLCN, QICN, CLCN
real, dimension(:,:,:), pointer :: ALH
real, dimension(: ), pointer :: PREF
real, dimension(IM,JM,1:LM-1) :: TVE, RDZ
real, dimension(IM,JM,LM) :: THV, TV, Z, DMI, PLO, QS, DQS, QL, QI, QA
real, dimension(IM,JM,0:LM) :: PKE, ZLE
real, dimension(:,:,:), pointer :: EKH, EKM, KHLS, KMLS, KHRAD, KHSFC
real, dimension(:,: ), pointer :: BSTAR, USTAR, ZPBL, PPBL, WERAD, WESFC,VSCRAD,KERAD,DBUOY,ZSML,ZCLD,ZRADML,FRLAND
real, dimension(:,: ), pointer :: WEBRV,VSCBRV,DSIEMS,CHIS,ZCLDTOP,DELSINV,SMIXT,ZRADBS,CLDRF,VSCSFC,RADRCODE
real, dimension(:,:,:), pointer :: AKSODT, CKSODT
real, dimension(:,:,:), pointer :: AKQODT, CKQODT
real, dimension(:,:,:), pointer :: AKVODT, CKVODT
logical, dimension(IM,JM ) :: CONVECT
real :: LOUIS
real :: LAMBDAM, LAMBDAM2
real :: LAMBDAH, LAMBDAH2
real :: ZKMENV, ZKHENV
real :: MINTHICK
real :: MINSHEAR
real :: AKHMMAX
real :: C_B, LAMBDA_B, ZMAX_B,LOUIS_MEMORY
real :: PRANDTLSFC,PRANDTLRAD,BETA_RAD,BETA_SURF,KHRADFAC,TPFAC_SURF,ENTRATE_SURF
real :: PCEFF_SURF, KHSFCFAC,ZCHOKE
integer :: L,LOCK_ON
logical :: KH_alloc
logical :: KM_alloc
logical :: RI_alloc
logical :: DU_alloc
logical :: EKH_alloc
logical :: EKM_alloc
logical :: KHSFC_alloc
logical :: KHRAD_alloc
logical :: LWCRT_alloc
logical :: ZRADML_alloc
logical :: ZRADBS_alloc
logical :: ZSML_alloc
logical :: ZCLD_alloc
call MAPL_TimerOn(MAPL,"---PRELIMS")
! Get Sounding from the import state
!-----------------------------------
call MAPL_GetPointer(IMPORT, T, 'T', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, Q, 'QV', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, TH, 'TH', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, U, 'U', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, V, 'V', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT,VARFLT, 'VARFLT', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, PREF, 'PREF', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, RADLW, 'RADLW', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT,RADLWC, 'RADLWC', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, QLLS, 'QLLS', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, QILS, 'QILS', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, QLCN, 'QLCN', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, QICN, 'QICN', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, CLLS, 'CLLS', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, CLCN, 'CLCN', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, BSTAR, 'BSTAR', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT, USTAR, 'USTAR', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(IMPORT,FRLAND, 'FRLAND', RC=STATUS); VERIFY_(STATUS)
! Compute the edge heights using Arakawa-Suarez hydrostatic equation
!---------------------------------------------------------------------------
PKE = (PLE/MAPL_P00)**MAPL_KAPPA
ZLE(:,:,LM) = 0.0
do L=LM,1,-1
ZLE(:,:,L-1) = ZLE(:,:,L) + &
(MAPL_CP/MAPL_GRAV)*TH(:,:,L)*(PKE(:,:,L)-PKE(:,:,L-1))
end do
! Layer height, pressure, and virtual temperatures
!-------------------------------------------------
! First add up clouds to get total fraction, ice, and liquid
QL = QLCN + QLLS
QI = QICN + QILS
QA = CLCN + CLLS
Z = 0.5*(ZLE(:,:,0:LM-1)+ZLE(:,:,1:LM))
PLO = 0.5*(PLE(:,:,0:LM-1)+PLE(:,:,1:LM))
TV = T *( 1.0 + MAPL_VIREPS *Q - QL - QI )
THV = TV*(TH/T)
TVE = (TV (:,:,1:LM-1) + TV (:,:,2:LM))*0.5
! Miscellaneous factors
!----------------------
RDZ = PLE(:,:,1:LM-1) / ( MAPL_RGAS * TVE )
RDZ = RDZ / (Z(:,:, 1:LM-1)-Z(:,:, 2:LM))
DMI = (MAPL_GRAV*DT)/(PLE(:,:,1:LM)-PLE(:,:,0:LM-1))
! Get turbulence parameters from configuration
!---------------------------------------------
call ESMF_ConfigGetAttribute(CF,LOUIS ,Label=trim(COMP_NAME)//"_LOUIS:" ,default=5.0 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,LAMBDAM ,Label=trim(COMP_NAME)//"_LAMBDAM:" ,default=160.0 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,LAMBDAM2 ,Label=trim(COMP_NAME)//"_LAMBDAM2:",default=160.0 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,LAMBDAH ,Label=trim(COMP_NAME)//"_LAMBDAH:" ,default=160.0 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,LAMBDAH2 ,Label=trim(COMP_NAME)//"_LAMBDAH2:",default=160.0 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,ZKMENV ,Label=trim(COMP_NAME)//"_ZKMENV:" ,default=3000. ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,ZKHENV ,Label=trim(COMP_NAME)//"_ZKHENV:" ,default=3000. ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,MINTHICK ,Label=trim(COMP_NAME)//"_MINTHICK:",default=0.1 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,MINSHEAR ,Label=trim(COMP_NAME)//"_MINSHEAR:",default=0.0030 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,C_B ,Label=trim(COMP_NAME)//"_C_B:" ,default=2.50E-6 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,LAMBDA_B ,Label=trim(COMP_NAME)//"_LAMBDA_B:",default=1500. ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,AKHMMAX ,Label=trim(COMP_NAME)//"_AKHMMAX:" ,default=100. ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,LOCK_ON ,Label=trim(COMP_NAME)//"_LOCK_ON:" ,default=1 ,RC=STATUS)
call ESMF_ConfigGetAttribute(CF,PRANDTLSFC , Label=trim(COMP_NAME)//"_PRANDTLSFC:", default=1.0 , RC=STATUS)
call ESMF_ConfigGetAttribute(CF,PRANDTLRAD , Label=trim(COMP_NAME)//"_PRANDTLRAD:", default=0.75, RC=STATUS)
call ESMF_ConfigGetAttribute(CF,BETA_RAD , Label=trim(COMP_NAME)//"_BETA_RAD:", default=0.23, RC=STATUS)
call ESMF_ConfigGetAttribute(CF,BETA_SURF , Label=trim(COMP_NAME)//"_BETA_SURF:", default=0.23, RC=STATUS)
call ESMF_ConfigGetAttribute(CF,KHRADFAC , Label=trim(COMP_NAME)//"_KHRADFAC:", default=0.85, RC=STATUS)
call ESMF_ConfigGetAttribute(CF,KHSFCFAC , Label=trim(COMP_NAME)//"_KHSFCFAC:", default=0.85, RC=STATUS)
call ESMF_ConfigGetAttribute(CF,TPFAC_SURF, Label=trim(COMP_NAME)//"_TPFAC_SURF:", default=10.0, RC=STATUS)
call ESMF_ConfigGetAttribute(CF,ENTRATE_SURF,Label=trim(COMP_NAME)//"_ENTRATE_SURF:", default=1.e-3, RC=STATUS)
call ESMF_ConfigGetAttribute(CF,PCEFF_SURF, Label=trim(COMP_NAME)//"_PCEFF_SURF:", default=0.05, RC=STATUS)
call ESMF_ConfigGetAttribute(CF,LOUIS_MEMORY,Label=trim(COMP_NAME)//"_LOUIS_MEMORY:", default=-999., RC=STATUS)
! Get pointers from export state...
!-----------------------------------
call MAPL_GetPointer(EXPORT, KH, 'KH', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KM, 'KM', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, RI, 'RI', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, DU, 'DU', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, EKH, 'EKH', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, EKM, 'EKM', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KHLS, 'KHLS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KMLS, 'KMLS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KHSFC, 'KHSFC', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KHRAD, 'KHRAD', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, ZPBL, 'ZPBL', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, PPBL, 'PPBL', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, LWCRT, 'LWCRT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, WERAD, 'WERAD', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, WESFC, 'WESFC', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, DBUOY, 'DBUOY', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, VSCRAD,'VSCRAD',RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, VSCsfc,'VSCSFC',RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KERAD, 'KERAD', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, VSCBRV,'VSCBRV',RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, WEBRV, 'WEBRV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, CHIS, 'CHIS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, DSIEMS,'DSIEMS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, ZCLD, 'ZCLD', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, ZSML, 'ZSML', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, ZRADML, 'ZRADML', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, ZRADBS, 'ZRADBS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, ZCLDTOP,'ZCLDTOP', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, DELSINV,'DELSINV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, RADRCODE,'RADRCODE', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, SMIXT ,'SMIXT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, CLDRF ,'CLDRF', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, ALH ,'ALH', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, AKSODT, 'AKSODT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, CKSODT, 'CKSODT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, AKQODT, 'AKQODT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, CKQODT, 'CKQODT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, AKVODT, 'AKVODT', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, CKVODT, 'CKVODT', RC=STATUS)
VERIFY_(STATUS)
! Check which ones have to be allocated
!--------------------------------------
KH_alloc = .not.associated(KH )
KM_alloc = .not.associated(KM )
RI_alloc = .not.associated(RI )
DU_alloc = .not.associated(DU )
EKH_alloc = .not.associated(EKH )
EKM_alloc = .not.associated(EKM )
KHSFC_alloc = .not.associated(KHSFC)
KHRAD_alloc = .not.associated(KHRAD)
LWCRT_alloc = .not.associated(LWCRT)
ZRADML_alloc = .not.associated(ZRADML)
ZRADBS_alloc = .not.associated(ZRADBS)
ZCLD_alloc = .not.associated(ZCLD)
ZSML_alloc = .not.associated(ZSML)
! Allocate those
!---------------
if( KH_alloc) allocate(KH (IM,JM,0:LM))
if( KM_alloc) allocate(KM (IM,JM,0:LM))
if( RI_alloc) allocate(RI (IM,JM,0:LM))
if( DU_alloc) allocate(DU (IM,JM,0:LM))
if( EKH_alloc) allocate(EKH (IM,JM,0:LM))
if( EKM_alloc) allocate(EKM (IM,JM,0:LM))
if(KHSFC_alloc) allocate(KHSFC(IM,JM,0:LM))
if(KHRAD_alloc) allocate(KHRAD(IM,JM,0:LM))
if(LWCRT_alloc) allocate(LWCRT(IM,JM,1:LM))
if(ZRADML_alloc) allocate(ZRADML(IM,JM))
if(ZRADBS_alloc) allocate(ZRADBS(IM,JM))
if(ZCLD_alloc) allocate(ZCLD(IM,JM))
if(ZSML_alloc) allocate(ZSML(IM,JM))
! Need DQSAT for Lock scheme
DQS = GEOS_DQSAT(T, PLO, qsat=QS, PASCALS=.true. )
!===> Running 1-2-1 smooth of bottom 5 levels of Virtual Pot. Temp.
THV(:,:,LM) = THV(:,:,LM-1)*0.25 + THV(:,:,LM)*0.75
do L = LM-1,LM-5,-1
THV(:,:,L) = THV(:,:,L-1)*0.25 + THV(:,:,L)*0.50 + THV(:,:,L+1)*0.25
end do
call MAPL_TimerOff(MAPL,"---PRELIMS")
! Refresh diffusivities: First compute Louis...
! ---------------------------------------------
call MAPL_TimerOn(MAPL,"---LOUIS")
call LOUIS_KS(KH(:,:,1:LM-1), KM(:,:,1:LM-1), &
RI(:,:,1:LM-1), DU(:,:,1:LM-1), &
Z, ZLE(:,:,1:LM-1), THV, U, V, &
LOUIS, MINSHEAR, MINTHICK, &
LAMBDAM, LAMBDAM2, LAMBDAH, LAMBDAH2, &
ZKMENV, ZKHENV, &
AKHMMAX, &
ALH, &
RC=STATUS )
VERIFY_(STATUS)
if(associated(KHLS)) KHLS = KH
if(associated(KMLS)) KMLS = KM
call MAPL_TimerOff(MAPL,"---LOUIS")
! ...then add Lock.
!--------------------
LWCRT = RADLW - RADLWC
if (LOCK_ON==1) then
call MAPL_TimerOn(MAPL,"---LOCK")
CALL ENTRAIN(1,IM,1,JM, RADLW,USTAR,BSTAR,FRLAND, &
T,Q,QL,QI,QA,U,V,Z,PLO,ZLE,PLE, &
QS,DQS, &
KM (:,:,0:LM-1), &
KH (:,:,0:LM-1), &
EKM (:,:,0:LM-1), &
EKH (:,:,0:LM-1), &
KHSFC(:,:,0:LM-1), &
KHRAD(:,:,0:LM-1), &
ZCLD, &
ZRADML, &
ZRADBS, &
ZSML, &
! these are optionally allocated
ZCLDTOP, &
WESFC, &
WERAD, &
DBUOY, &
VSCSFC, &
VSCRAD, &
KERAD, &
VSCBRV, &
WEBRV, &
DSIEMS, &
CHIS, &
DELSINV, &
SMIXT, &
CLDRF, &
RADRCODE, &
! these are input params
PRANDTLSFC,PRANDTLRAD, &
BETA_SURF,BETA_RAD, &
TPFAC_SURF,ENTRATE_SURF,PCEFF_SURF, &
KHRADFAC, KHSFCFAC )
call MAPL_TimerOff(MAPL,"---LOCK")
else
EKM (:,:,0:LM-1) = 0.0
EKH (:,:,0:LM-1) = 0.0
KHSFC(:,:,0:LM-1) = 0.0
KHRAD(:,:,0:LM-1) = 0.0
endif
! PBL-top diagnostic
! -----------------------------------------
if (associated(ZPBL)) then
ZPBL = MAPL_UNDEF
where( (KH(:,:,LM-1) < 2.) ) ZPBL = Z(:,:,LM)
do L=LM-2,1,-1
where( (KH(:,:,L) < 2.) .and. (KH(:,:,L+1) >= 2.) .and. (ZPBL == MAPL_UNDEF ) )
ZPBL = Z(:,:,L+1)
endwhere
enddo
endif
if (associated(PPBL)) then
PPBL = MAPL_UNDEF
where( (KH(:,:,LM-1) < 2.) ) PPBL = PLO(:,:,LM)
do L=LM-2,1,-1
where( (KH(:,:,L) < 2.) .and. (KH(:,:,L+1) >= 2.) .and. (PPBL == MAPL_UNDEF ) )
PPBL = PLO(:,:,L+1)
endwhere
enddo
endif
! Second difference coefficients for scalars; RDZ is RHO/DZ, DMI is (G DT)/DP
! ---------------------------------------------------------------------------
CKS(:,:,1:LM-1) = -KH(:,:,1:LM-1) * RDZ(:,:,1:LM-1)
AKS(:,:,1 ) = 0.0
AKS(:,:,2:LM ) = CKS(:,:,1:LM-1) * DMI(:,:,2:LM )
CKS(:,:,1:LM-1) = CKS(:,:,1:LM-1) * DMI(:,:,1:LM-1)
CKS(:,:, LM ) = - CT * DMI(:,:, LM )
! Water vapor can differ at the surface
!--------------------------------------
CKQ = CKS
AKQ = AKS
CKQ(:,:, LM ) = - CQ * DMI(:,:, LM )
! Second difference coefficients for winds
! EKV is saved to use in the frictional heating calc.
! ---------------------------------------------------
EKV(:,:,1:LM-1) = -KM(:,:,1:LM-1) * RDZ(:,:,1:LM-1)
AKV(:,:,1 ) = 0.0
AKV(:,:,2:LM ) = EKV(:,:,1:LM-1) * DMI(:,:,2:LM )
CKV(:,:,1:LM-1) = EKV(:,:,1:LM-1) * DMI(:,:,1:LM-1)
EKV(:,:,1:LM-1) = -MAPL_GRAV * EKV(:,:,1:LM-1)
CKV(:,:, LM ) = - CU * DMI(:,:, LM )
EKV(:,:, LM ) = MAPL_GRAV * CU
! Setup the tridiagonal matrix
! ----------------------------
BKS = 1.00 - (AKS+CKS)
BKQ = 1.00 - (AKQ+CKQ)
BKV = 1.00 - (AKV+CKV)
! Add the topographic roughness term
!-----------------------------------
if( associated(AKSODT) ) then
AKSODT = -AKS/DT
AKSODT(:,:, 1) = 0.0
endif
if( associated(CKSODT) ) then
CKSODT = -CKS/DT
CKSODT(:,:,LM) = 0.0
endif
if( associated(AKQODT) ) then
AKQODT = -AKQ/DT
AKQODT(:,:, 1) = 0.0
endif
if( associated(CKQODT) ) then
CKQODT = -CKQ/DT
CKQODT(:,:,LM) = 0.0
endif
if( associated(AKVODT) ) AKVODT = -AKV/DT
if( associated(CKVODT) ) CKVODT = -CKV/DT
!BOP
!
! Orograpghic drag follows Beljaars (2003):
! $$
! \frac{\partial}{\partial z}\frac{\tau}{\rho} = \frac{C_B}{\lambda_B} |U(z)| U(z)
! e^{-\tilde{z}^\frac{3}{2}}\tilde{z}^{-1.2},
! $$
! where $z$ is the height above the surface in meters,
! $\tilde{z}=\frac{z}{\lambda_B}$, $\tau$ is the orographic stress at $z$,
! $\rho$ is the air density, $U(z)$ is the wind velocity, and $\lambda_B$ is a vertical length scale.
! Beljaars uses $\lambda_B = 1500$m, for which the non-dimensional parameter $C_B = 2.50 \times 10^{-6}$.
! These are the default values, but both can be modified from the configuration. To avoid underflow.
! the tendency is set to zero once $\tilde{z}$ exceeds 4 (i.e., 6 km from the surface for default values).
!
!EOP
call MAPL_TimerOn(MAPL,"---BELJAARS")
do L=LM,1,-1
where(Z(:,:,L) < 4.0*LAMBDA_B)
FKV(:,:,L) = Z(:,:,L)*(1.0/LAMBDA_B)
FKV(:,:,L) = VARFLT * exp(-FKV(:,:,L)*sqrt(FKV(:,:,L)))*(FKV(:,:,L)**(-1.2))
FKV(:,:,L) = (C_B/LAMBDA_B)*sqrt(U(:,:,L)**2+V(:,:,L)**2)*FKV(:,:,L)
elsewhere
FKV(:,:,L) = 0.0
end where
BKV(:,:,L) = BKV(:,:,L) + DT*FKV(:,:,L)
FKV(:,:,L) = FKV(:,:,L) * (PLE(:,:,L)-PLE(:,:,L-1))
end do
call MAPL_TimerOff(MAPL,"---BELJAARS")
call MAPL_TimerOn(MAPL,"---DECOMP")
! Do LU decomposition; C is not modified.
! On exit, B is the main diagonals of the LU
! decomposition, and A is the r.h.s multiplier.
!----------------------------------------------
call VTRILU(AKS,BKS,CKS)
call VTRILU(AKQ,BKQ,CKQ)
call VTRILU(AKV,BKV,CKV)
! Get the sensitivity of solution to a unit
! change in the surface value. B and C are
! not modified.
!------------------------------------------
call VTRISOLVESURF(BKS,CKS,DKS)
call VTRISOLVESURF(BKQ,CKQ,DKQ)
call VTRISOLVESURF(BKV,CKV,DKV)
call MAPL_TimerOff(MAPL,"---DECOMP")
! Deallocate the temporaries
!---------------------------
if(KH_alloc ) deallocate(KH)
if(KM_alloc ) deallocate(KM)
if(RI_alloc ) deallocate(RI)
if(DU_alloc ) deallocate(DU)
if(EKH_alloc ) deallocate(EKH)
if(EKM_alloc ) deallocate(EKM)
if(KHSFC_alloc) deallocate(KHSFC)
if(KHRAD_alloc) deallocate(KHRAD)
if(LWCRT_alloc) deallocate(LWCRT)
if(ZRADML_alloc) deallocate(ZRADML)
if(ZRADBS_alloc) deallocate(ZRADBS)
if(ZCLD_alloc) deallocate(ZCLD)
if(ZSML_alloc) deallocate(ZSML)
RETURN_(ESMF_SUCCESS)
end subroutine REFRESH
!=============================================================================
!=============================================================================
!BOP
! !CROUTINE: DIFFUSE -- Solves for semi-implicit diffusive tendencies assuming fixed surface conditions.
! !INTERFACE:
subroutine DIFFUSE(IM,JM,LM,RC) 1,6
! !ARGUMENTS:
integer, intent(IN) :: IM,JM,LM
integer, optional, intent(OUT) :: RC
! !DESCRIPTION: {\tt DIFFUSE} computes semi-implicit tendencies of all fields in
! the TR bundle. Each field is examined for three attributes: {\tt DiffuseLike},
! {\tt FriendlyToTURBULENCE}, and {\tt WeightedTendency}. These determine the behavior of
! {\tt DIFFUSE} for that field. {\tt DiffuseLike} can be either 'U', 'Q', or 'S'; the default is 'Q'.
! {\tt FriendlyToTURBULENCE}, and {\tt WeightedTendency} are ESMF logicals.
! If {\tt FriendlyToTURBULENCE} is true, the field in TR is updated directly; otherwise
! it is left untouched. In either case, If the corresponding pointer TRI bundle is associated, the
! tendencies are returned there. If {\tt WeightedTendency} is true, the tendency in TRI, if any,
! is pressure weighted.
!EOP
character(len=ESMF_MAXSTR) :: IAm='Diffuse'
integer :: STATUS
character(len=ESMF_MAXSTR) :: TYPE
type (ESMF_Field) :: FIELD
type (ESMF_Array) :: ARRAY
type (ESMF_Bundle) :: TR
type (ESMF_Bundle) :: TRI
type (ESMF_Bundle) :: TRG
type (ESMF_Bundle) :: FSTAR
type (ESMF_Bundle) :: DFSTAR
real, dimension(:,:,:), pointer :: S, SOI, SOD
real, dimension(:,:), pointer :: SG, SF, SDF, CX, SRG
real, dimension(:,:,:), pointer :: DX
real, dimension(:,:,:), pointer :: AK, BK, CK
integer :: KM, K,L
type(ESMF_logical) :: FRIENDLY
type(ESMF_logical) :: WEIGHTED
real, dimension(IM,JM,LM) :: DP, SX
! Get the bundles containing the quantities to be diffused,
! their tendencies, their surface values, their surface
! fluxes, and the derivatives of their surface fluxes
! wrt the surface values.
!----------------------------------------------------------
call ESMF_StateGetBundle(IMPORT, 'TR' , TR, RC=STATUS); VERIFY_(STATUS)
call ESMF_StateGetBundle(IMPORT, 'TRG', TRG, RC=STATUS); VERIFY_(STATUS)
call ESMF_StateGetBundle(EXPORT, 'TRI', TRI, RC=STATUS); VERIFY_(STATUS)
call ESMF_StateGetBundle(EXPORT, 'FSTAR', FSTAR, RC=STATUS); VERIFY_(STATUS)
call ESMF_StateGetBundle(EXPORT, 'DFSTAR', DFSTAR, RC=STATUS); VERIFY_(STATUS)
! Count the firlds in TR...
!--------------------------
call ESMF_BundleGet(TR, fieldCOUNT=KM, RC=STATUS)
VERIFY_(STATUS)
! ...and make sure the other bundles are the same.
!-------------------------------------------------
call ESMF_BundleGet(TRI, FieldCount=K , RC=STATUS)
VERIFY_(STATUS)
ASSERT_(KM==K)
call ESMF_BundleGet(TRG, FieldCount=K , RC=STATUS)
VERIFY_(STATUS)
ASSERT_(KM==K)
call ESMF_BundleGet(FSTAR, FieldCount=K , RC=STATUS)
VERIFY_(STATUS)
ASSERT_(KM==K)
call ESMF_BundleGet(DFSTAR, FieldCount=K , RC=STATUS)
VERIFY_(STATUS)
ASSERT_(KM==K)
! Pressure thickness of layers
!-----------------------------
DP = PLE(:,:,1:LM)-PLE(:,:,0:LM-1)
! Loop over all quantities to be diffused.
!----------------------------------------
do K=1,KM
! Get the Kth Field from tracer bundle
!-------------------------------------
call ESMF_BundleGetField (TR, K, FIELD, RC=STATUS)
VERIFY_(STATUS)
! Get item's diffusion type (U, S or Q; default is Q)
!----------------------------------------------------
call ESMF_FieldGetAttribute (FIELD, NAME="DiffuseLike", VALUE=TYPE, RC=STATUS)
if(STATUS /= ESMF_SUCCESS) TYPE = 'Q'
! Get item's friendly status (default is not friendly)
!-----------------------------------------------------
call ESMF_FieldGetAttribute (FIELD, NAME="FriendlyToTURBULENCE",VALUE=FRIENDLY, RC=STATUS)
if(STATUS /= ESMF_SUCCESS) FRIENDLY = ESMF_false
! Get item's weighting (default is unweighted tendencies)
!--------------------------------------------------------
call ESMF_FieldGetAttribute (FIELD, NAME="WeightedTendency", VALUE=WEIGHTED, RC=STATUS)
if(STATUS /= ESMF_SUCCESS) WEIGHTED = ESMF_false
! Get pointer to the quantity, its tendency, its surface value,
! the surface flux, and the sensitivity of the surface flux.
! -------------------------------------------------------------
call ESMFL_BundleGetPointerToData(TR , K, S , RC=STATUS); VERIFY_(STATUS)
call ESMFL_BundleGetPointerToData(TRI , K, SOI, RC=STATUS); VERIFY_(STATUS)
call ESMFL_BundleGetPointerToData(TRG , K, SRG, RC=STATUS); VERIFY_(STATUS)
call ESMFL_BundleGetPointerToData(FSTAR , K, SF , RC=STATUS); VERIFY_(STATUS)
call ESMFL_BundleGetPointerToData(DFSTAR, K, SDF, RC=STATUS); VERIFY_(STATUS)
! The quantity must exist; others are optional.
!----------------------------------------------
ASSERT_(associated(S ))
! If the surface values does not exists, we assume zero flux.
!------------------------------------------------------------
if(associated(SRG)) then
SG => SRG
else
allocate (SG(0,0), stat=STATUS)
VERIFY_(STATUS)
end if
! Pick the right exchange coefficients
!-------------------------------------
if ( TYPE=='U' ) then ! Momentum
CX => CU
DX => DKV
AK => AKV; BK => BKV; CK => CKV
else if( TYPE=='Q' ) then ! Water Vapor
CX => CQ
DX => DKQ
AK => AKQ; BK => BKQ; CK => CKQ
else if( TYPE=='S' ) then ! Heat or other scalars
CX => CT
DX => DKS
AK => AKS; BK => BKS; CK => CKS
else
RETURN_(ESMF_FAILURE)
endif
! Copy diffused quantity to temp buffer
! ------------------------------------------
SX = S
! Solve for semi-implicit changes. This modifies SX
! -------------------------------------------------
call VTRISOLVE(AK,BK,CK,SX,SG)
! Compute the surface fluxes
!---------------------------
if(associated(SF)) then
if(size(SG)>0) then
SF = CX*(SG - SX(:,:,LM))
else
SF = 0.0
end if
end if
! Create tendencies
!------------------
if(associated(SOI)) then
if( WEIGHTED==ESMF_TRUE ) then
SOI = ( (SX - S)/DT )*DP
else
SOI = ( (SX - S)/DT )
endif
end if
! Update friendlies
!------------------
if(FRIENDLY==ESMF_TRUE) then
S = SX
end if
! Compute the derivative of the surface flux wrt the surface value
!-----------------------------------------------------------------
if(associated(SDF)) then
SDF = CX * (1.0-DX(:,:,LM))
endif
if(.not.associated(SRG)) then
deallocate (SG)
end if
enddo ! End loop over all quantities to be diffused
! -----------------------------------------------------
RETURN_(ESMF_SUCCESS)
end subroutine DIFFUSE
end subroutine RUN1
!*********************************************************************
!*********************************************************************
!*********************************************************************
!BOP
! !IROUTINE: RUN2 -- The second run stage for the TURBULENCE component
! !INTERFACE:
subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ),179
! !ARGUMENTS:
type(ESMF_GridComp), intent(inout) :: GC ! Gridded component
type(ESMF_State), intent(inout) :: IMPORT ! Import state
type(ESMF_State), intent(inout) :: EXPORT ! Export state
type(ESMF_Clock), intent(inout) :: CLOCK ! The clock
integer, optional, intent( out) :: RC ! Error code:
! !DESCRIPTION: Second run stage of {\tt GEOS\_TurbulenceGridComp} performs
! the updates due to changes in surface quantities. Its input are the changes in
! surface quantities during the time step. It can also compute the frictional
! dissipation terms as exports, but these are not added to the temperatures.
!EOP
! ErrLog Variables
character(len=ESMF_MAXSTR) :: IAm
integer :: STATUS
character(len=ESMF_MAXSTR) :: COMP_NAME
! Local derived type aliases
type (MAPL_MetaComp), pointer :: MAPL
type (ESMF_Config ) :: CF
type (ESMF_State ) :: INTERNAL
! Local variables
integer :: IM, JM, LM
real :: DT
! Begin...
!---------
! Get my name and set-up traceback handle
! ---------------------------------------
call ESMF_GridCompGet( GC, NAME=COMP_NAME, RC=STATUS )
VERIFY_(STATUS)
Iam = trim(COMP_NAME) // 'Run2'
! Get my internal MAPL_Generic state
!-----------------------------------
call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerOn(MAPL,"TOTAL")
call MAPL_TimerOn(MAPL,"-RUN2")
! Get parameters from generic state.
!-----------------------------------
call MAPL_Get(MAPL, &
IM=IM, JM=JM, LM=LM, &
INTERNAL_ESMF_STATE=INTERNAL, &
RC=STATUS )
VERIFY_(STATUS)
! Get configuration from component
!---------------------------------
call ESMF_GridCompGet( GC, CONFIG = CF, RC=STATUS )
VERIFY_(STATUS)
! Get application's timestep from configuration
!----------------------------------------------
call ESMF_ConfigGetAttribute( CF, DT, Label="RUN_DT:" , RC=STATUS)
VERIFY_(STATUS)
! Solve the free atmosphere problem
! ---------------------------------
call MAPL_TimerOn (MAPL,"--UPDATE")
call UPDATE(IM,JM,LM,RC=STATUS)
VERIFY_(STATUS)
call MAPL_TimerOff(MAPL,"--UPDATE")
! All done with RUN
!-------------------
call MAPL_TimerOff(MAPL,"-RUN2")
call MAPL_TimerOff(MAPL,"TOTAL")
RETURN_(ESMF_SUCCESS)
contains
!BOP
! !CROUTINE: UPDATE -- Updates diffusive effects for changes at surface.
! !INTERFACE:
subroutine UPDATE(IM,JM,LM,RC) 8,13
! !ARGUMENTS:
integer, intent(IN) :: IM,JM,LM
integer, optional, intent(OUT) :: RC
! !DESCRIPTION:
! Some description
!EOP
character(len=ESMF_MAXSTR) :: IAm='Update'
integer :: STATUS
character(len=ESMF_MAXSTR) :: TYPE
type (ESMF_Field) :: FIELD
type (ESMF_Bundle) :: TR
type (ESMF_Bundle) :: TRI
type (ESMF_Bundle) :: DTG
type (ESMF_Bundle) :: FSTAR
type (ESMF_Bundle) :: DFSTAR
real, dimension(:,:,:), pointer :: PLE
real, dimension(:,:,:), pointer :: S, SOI, INTDIS, TOPDIS
real, dimension(:,: ), pointer :: DSG, SF, SDF, SRFDIS
real, dimension(:,: ), pointer :: KETRB, KESRF, KETOP, KEINT
real, dimension(:,:,:), pointer :: DKS, DKV, DKQ, DKX, EKV, FKV
integer :: KM, K, L
type(ESMF_logical) :: FRIENDLY
type(ESMF_logical) :: WEIGHTED
real, dimension(IM,JM,LM) :: DP, SX
real, dimension(IM,JM,LM-1) :: DF
! Pressure-weighted dissipation heating rates
!--------------------------------------------
call MAPL_GetPointer(EXPORT, KETRB , 'KETRB' , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KESRF , 'KESRF' , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KETOP , 'KETOP' , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, KEINT , 'KEINT' , RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, SRFDIS, 'SRFDIS', &
alloc=associated(KETRB) .or. associated(KESRF), &
RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, INTDIS, 'INTDIS', &
alloc=associated(KETRB) .or. associated(KEINT), &
RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, TOPDIS, 'TOPDIS', &
alloc=associated(KETRB) .or. associated(KETOP), &
RC=STATUS)
VERIFY_(STATUS)
! Get imports
!------------
call MAPL_GetPointer(IMPORT, PLE, 'PLE', RC=STATUS); VERIFY_(STATUS)
! Get the tendecy sensitivities computed in RUN1
!-----------------------------------------------
call MAPL_GetPointer(INTERNAL, DKS, 'DKS', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, DKV, 'DKV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, DKQ, 'DKQ', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, EKV, 'EKV', RC=STATUS)
VERIFY_(STATUS)
call MAPL_GetPointer(INTERNAL, FKV, 'FKV', RC=STATUS)
VERIFY_(STATUS)
! Get the bundles containing the quantities to be diffused,
! their tendencies, their surface values, their surface
! fluxes, and the derivatives of their surface fluxes
! wrt the surface values.
!----------------------------------------------------------
call ESMF_StateGetBundle(IMPORT, 'TR' , TR, RC=STATUS); VERIFY_(STATUS)
call ESMF_StateGetBundle(IMPORT, 'DTG', DTG, RC=STATUS); VERIFY_(STATUS)
call ESMF_StateGetBundle(EXPORT, 'TRI', TRI, RC=STATUS); VERIFY_(STATUS)
call ESMF_StateGetBundle(EXPORT, 'FSTAR' , FSTAR, RC=STATUS); VERIFY_(STATUS)
call ESMF_StateGetBundle(EXPORT, 'DFSTAR', DFSTAR, RC=STATUS); VERIFY_(STATUS)
! Count them...
!--------------
call ESMF_BundleGet(TR , FieldCount=KM, RC=STATUS)
VERIFY_(STATUS)
! and make sure the other bundles are the same.
!----------------------------------------------
call ESMF_BundleGet(DTG, FieldCount=K , RC=STATUS)
VERIFY_(STATUS)
ASSERT_(KM==K)
! Clear the accumulators for the dissipation.
!--------------------------------------------
if(associated(SRFDIS)) SRFDIS = 0.0
if(associated(INTDIS)) INTDIS = 0.0
if(associated(TOPDIS)) TOPDIS = 0.0
if(associated(KETRB )) KETRB = 0.0
if(associated(KESRF )) KESRF = 0.0
if(associated(KETOP )) KETOP = 0.0
if(associated(KEINT )) KEINT = 0.0
! Pressure thickness of layers
!-----------------------------
DP = PLE(:,:,1:LM)-PLE(:,:,0:LM-1)
! Loop over all quantities to be diffused.
!-----------------------------------------
do K=1,KM
! Get Kth field from bundle
!--------------------------
call ESMF_BundleGetField (TR, K, FIELD, RC=STATUS)
VERIFY_(STATUS)
! Get item's diffusion type (U, S or Q; default is Q)
!----------------------------------------------------
call ESMF_FieldGetAttribute (FIELD, NAME="DiffuseLike", VALUE=TYPE, RC=STATUS)
if(STATUS /= ESMF_SUCCESS) TYPE = 'Q'
! Get item's friendly status (default is FALSE)
!----------------------------------------------
call ESMF_FieldGetAttribute (FIELD, NAME="FriendlyToTURBULENCE",VALUE=FRIENDLY, RC=STATUS)
if(STATUS /= ESMF_SUCCESS) FRIENDLY = ESMF_false
! Get item's weighting attribute (default is TRUE)
!-------------------------------------------------
call ESMF_FieldGetAttribute (FIELD, NAME="WeightedTendency", VALUE=WEIGHTED, RC=STATUS)
if(STATUS /= ESMF_SUCCESS) WEIGHTED = ESMF_false
! Get pointers to the quantity, its tendency, its surface increment,
! the preliminary surface flux, and the sensitivity of the surface
! flux to the surface value.
! ------------------------------------------------------------------
call ESMFL_BundleGetPointerToData(TR , K, S , RC=STATUS)
VERIFY_(STATUS)
call ESMFL_BundleGetPointerToData(TRI , K, SOI, RC=STATUS)
VERIFY_(STATUS)
call ESMFL_BundleGetPointerToData(DTG , K, DSG, RC=STATUS)
VERIFY_(STATUS)
call ESMFL_BundleGetPointerToData(FSTAR , K, SF , RC=STATUS)
VERIFY_(STATUS)
call ESMFL_BundleGetPointerToData(DFSTAR, K, SDF, RC=STATUS)
VERIFY_(STATUS)
! Point to the appropriate sensitivity
!--------------------------------------
if ( TYPE=='U' ) then
DKX => DKV
else if ( TYPE=='Q' ) then
DKX => DKQ
else if ( TYPE=='S' ) then
DKX => DKS
else
RETURN_(ESMF_FAILURE)
end if
! Update diffused quantity
!-------------------------
SX = S
if(associated(DSG)) then
do L=1,LM
SX(:,:,L) = SX(:,:,L) + DKX(:,:,L)*DSG
end do
end if
! Replace friendly
!-----------------
if(FRIENDLY==ESMF_TRUE) then
S = SX
end if
! Increment the dissipation
!--------------------------
if( TYPE=='U' ) then
if(associated(KETRB )) KETRB = 0.0
if(associated(KESRF )) KESRF = 0.0
if(associated(KETOP )) KETOP = 0.0
if(associated(KEINT )) KEINT = 0.0
if(associated(INTDIS)) then
DF = (0.5/(MAPL_CP))*EKV(:,:,1:LM-1)*(SX(:,:,1:LM-1)-SX(:,:,2:LM))**2
INTDIS(:,:,1:LM-1) = INTDIS(:,:,1:LM-1) + DF
INTDIS(:,:,2:LM ) = INTDIS(:,:,2:LM ) + DF
if(associated(KETRB)) then
do L=1,LM
KETRB = KETRB - INTDIS(:,:,L)* (MAPL_CP/MAPL_GRAV)
end do
end if
if(associated(KEINT)) then
do L=1,LM
KEINT = KEINT - INTDIS(:,:,L)* (MAPL_CP/MAPL_GRAV)
end do
end if
endif
if(associated(TOPDIS)) then
TOPDIS = TOPDIS + (1.0/(MAPL_CP))*FKV*SX**2
if(associated(KETRB)) then
do L=1,LM
KETRB = KETRB - TOPDIS(:,:,L)* (MAPL_CP/MAPL_GRAV)
end do
end if
if(associated(KETOP)) then
do L=1,LM
KETOP = KETOP - TOPDIS(:,:,L)* (MAPL_CP/MAPL_GRAV)
end do
end if
endif
if(associated(SRFDIS)) then
SRFDIS = SRFDIS + (1.0/(MAPL_CP))*EKV(:,:,LM)*SX(:,:,LM)**2
if(associated(KETRB)) KETRB = KETRB - SRFDIS* (MAPL_CP/MAPL_GRAV)
if(associated(KESRF)) KESRF = KESRF - SRFDIS* (MAPL_CP/MAPL_GRAV)
endif
end if
! Update tendencies
! -----------------
if(associated(SOI).and.associated(DSG)) then
if( WEIGHTED==ESMF_TRUE ) then
do L=1,LM
SOI(:,:,L) = SOI(:,:,L) + (DKX(:,:,L)*DSG/DT)*DP(:,:,L)
end do
else
do L=1,LM
SOI(:,:,L) = SOI(:,:,L) + (DKX(:,:,L)*DSG/DT)
end do
endif
end if
! Update surface fluxes
! ---------------------
if(associated(SF).and.associated(DSG)) then
SF = SF + DSG*SDF
end if
enddo ! End loop over all quantities to be diffused
!--------------------------------------------------------
RETURN_(ESMF_SUCCESS)
end subroutine UPDATE
end subroutine RUN2
!*********************************************************************
!*********************************************************************
!*********************************************************************
!BOP
! !IROUTINE: VTRILU -- Does LU decomposition of tridiagonal matrix.
! !INTERFACE:
subroutine VTRILU ( A,B,C ) 3
! !ARGUMENTS:
real, dimension(:,:,:), intent(IN ) :: C
real, dimension(:,:,:), intent(INOUT) :: A, B
! !DESCRIPTION: {\tt VTRILU} performs an $LU$ decomposition on
! a tridiagonal matrix $M=LU$.
!
! $$
! M = \left( \begin{array}{ccccccc}
! b_1 & c_1 & & & & & \\
! a_2 & b_2 & c_2 & & & & \\
! & \cdot& \cdot & \cdot & & & \\
! & & \cdot& \cdot & \cdot & & \\
! && & \cdot& \cdot & \cdot & \\
! &&&& a_{K-1} & b_{K-1} & c_{K-1} \\
! &&&&& a_{K} & b_{K}
! \end{array} \right)
! $$
!
! $$
! \begin{array}{lr}
! L = \left( \begin{array}{ccccccc}
! 1 &&&&&& \\
! \hat{a}_2 & 1 & &&&& \\
! & \cdot& \cdot & & & & \\
! & & \cdot& \cdot & && \\
! && & \cdot& \cdot & & \\
! &&&& \hat{a}_{K-1} & 1 & \\
! &&&&& \hat{a}_{K} & 1
! \end{array} \right)
! &
! U = \left( \begin{array}{ccccccc}
! \hat{b}_1 & c_1 &&&&& \\
! & \hat{b}_2 & c_2 &&&& \\
! & & \cdot & \cdot & & & \\
! & & & \cdot & \cdot && \\
! && & & \cdot & \cdot & \\
! &&&& & \hat{b}_{K-1} & c_{K-1} \\
! &&&&& & \hat{b}_{K}
! \end{array} \right)
! \end{array}
! $$
!
! On input, A, B, and C contain, $a_k$, $b_k$, and $c_k$
! the lower, main, and upper diagonals of the matrix, respectively.
! On output, B contains $1/\hat{b}_k$, the inverse of the main diagonal of $U$,
! and A contains $\hat{a}_k$,
! the lower diagonal of $L$. C contains the upper diagonal of the original matrix and of $U$.
!
! The new diagonals $\hat{a}_k$ and $\hat{b}_k$ are:
! $$
! \begin{array}{rcl}
! \hat{b}_1 & = & b_1, \\
! \hat{a}_k & = & \makebox[2 in][l]{$a_k / \hat{b}_{k-1}$,} k=2, K, \\
! \hat{b}_k & = & \makebox[2 in][l]{$b_k - c_{k-1} \hat{a}_k$,} k=2, K.
! \end{array}
! $$
!EOP
integer :: L
B(:,:,1) = 1. / B(:,:,1)
do L = 2,size(A,3)
A(:,:,L) = A(:,:,L) * B(:,:,L-1)
B(:,:,L) = 1. / ( B(:,:,L) - C(:,:,L-1) * A(:,:,L) )
enddo
return
end subroutine VTRILU
!*********************************************************************
!BOP
! !IROUTINE: VTRISOLVE -- Solves for tridiagonal system that has been decomposed by VTRILU
! !INTERFACE:
subroutine VTRISOLVE ( A,B,C,Y,YG ) 1
! !ARGUMENTS:
real, dimension(:,:,:), intent(IN ) :: A, B, C
real, dimension(:,:,:), intent(INOUT) :: Y
real, dimension(:,:), intent(IN) :: YG
! !DESCRIPTION: Solves tridiagonal system that has been LU decomposed
! $LU x = f$. This is done by first solving $L g = f$ for $g$, and
! then solving $U x = g$ for $x$. The solutions are:
! $$
! \begin{array}{rcl}
! g_1 & = & f_1, \\
! g_k & = & \makebox[2 in][l]{$f_k - g_{k-1} \hat{a}_{k}$,} k=2, K, \\
! \end{array}
! $$
! and
! $$
! \begin{array}{rcl}
! x_K & = & g_K /\hat{b}_K, \\
! x_k & = & \makebox[2 in][l]{($g_k - c_k g_{k+1}) / \hat{b}_{k}$,} k=K-1, 1 \\
! \end{array}
! $$
!
! On input A contains the $\hat{a}_k$, the lower diagonal of $L$,
! B contains the $1/\hat{b}_k$, inverse of the main diagonal of $U$,
! C contains the $c_k$, the upper diagonal of $U$. The forcing, $f_k$ is
!
! It returns the
! solution in the r.h.s input vector, Y. A has the multiplier from the
! decomposition, B the
! matrix (U), and C the upper diagonal of the original matrix and of U.
! YG is the LM+1 (Ground) value of Y.
!EOP
integer :: LM, L
LM = size(Y,3)
! Sweep down, modifying rhs with multiplier A
do L = 2,LM
Y(:,:,L) = Y(:,:,L) - Y(:,:,L-1) * A(:,:,L)
enddo
! Sweep up, solving for updated value. Note B has the inverse of the main diagonal
if(size(YG)>0) then
Y(:,:,LM) = (Y(:,:,LM) - C(:,:,LM) * YG )*B(:,:,LM)
else
Y(:,:,LM) = Y(:,:,LM)*B(:,:,LM-1)/(B(:,:,LM-1) - A(:,:,LM)*(1.0+C(:,:,LM-1)*B(:,:,LM-1) ))
endif
do L = LM-1,1,-1
Y(:,:,L) = (Y(:,:,L ) - C(:,:,L ) * Y(:,:,L+1))*B(:,:,L )
enddo
return
end subroutine VTRISOLVE
!*********************************************************************
!BOP
! !IROUTINE: VTRISOLVESURF -- Solves for sensitivity to surface value
! !INTERFACE:
subroutine VTRISOLVESURF( B,C,Y ) 3
! !ARGUMENTS:
real, dimension(:,:,:), intent(IN ) :: B, C
real, dimension(:,:,:), intent( OUT) :: Y
! !DESCRIPTION: Solves tridiagonal system that has been LU decomposed
! for the special case
! where the surface Y (YG) is 1 and the rest of the input Ys are 0.
! Everything else is as in {\tt VTRISOLVE}. This gives the sensitivity of the
! solution to a unit change in the surface values.
!EOP
integer :: LM, L
LM = size(B,3)
Y(:,:,LM) = -C(:,:,LM) * B(:,:,LM)
do L = LM-1,1,-1
Y(:,:,L) = -C(:,:,L) * Y(:,:,L+1) * B(:,:,L)
enddo
return
end subroutine VTRISOLVESURF
!*********************************************************************
!BOP
! !IROUTINE: GETKS -- Computes atmospheric diffusivities at interior levels
! !INTERFACE:
subroutine LOUIS_KS( & 1
KH,KM,RI,DU, &
ZZ,ZE,PV,UU,VV, &
LOUIS, MINSHEAR, MINTHICK, &
LAMBDAM, LAMBDAM2,LAMBDAH, LAMBDAH2, &
ZKMENV, ZKHENV, &
KHMMAX, &
ALH_DIAG, &
RC )
! !ARGUMENTS:
real, intent(IN ) :: ZZ(:,:,:) ! Height of layer center above the surface (m).
real, intent(IN ) :: PV(:,:,:) ! Virtual potential temperature at layer center (K).
real, intent(IN ) :: UU(:,:,:) ! Eastward velocity at layer center (m s-1).
real, intent(IN ) :: VV(:,:,:) ! Northward velocity at layer center (m s-1).
real, intent(IN ) :: ZE(:,:,:) ! Height of layer base above the surface (m).
real, intent(IN ) :: LOUIS ! Louis scheme parameters (usually 5).
real, intent(IN ) :: MINSHEAR ! Min shear allowed in Ri calculation (s-1).
real, intent(IN ) :: MINTHICK ! Min layer thickness (m).
real, intent(IN ) :: KHMMAX ! Maximum allowe diffusivity (m+2 s-1).
real, intent(IN ) :: LAMBDAM ! Blackadar(1962) length scale parameter for momentum (m).
real, intent(IN ) :: LAMBDAM2 ! Second Blackadar parameter for momentum (m).
real, intent(IN ) :: LAMBDAH ! Blackadar(1962) length scale parameter for heat (m).
real, intent(IN ) :: LAMBDAH2 ! Second Blackadar parameter for heat (m).
real, intent(IN ) :: ZKMENV ! Transition height for Blackadar param for momentum (m)
real, intent(IN ) :: ZKHENV ! Transition height for Blackadar param for heat (m)
real, intent(OUT) :: KM(:,:,:) ! Momentum diffusivity at base of each layer (m+2 s-1).
real, intent(OUT) :: KH(:,:,:) ! Heat diffusivity at base of each layer (m+2 s-1).
real, intent(OUT) :: RI(:,:,:) ! Richardson number
real, intent(OUT) :: DU(:,:,:) ! Magnitude of wind shear (s-1).
integer, optional, intent(OUT) :: RC
real, pointer :: ALH_DIAG(:,:,:) ! Blackadar Length Scale diagnostic (m) [Optional]
! !DESCRIPTION: Computes Louis et al.(1979) Richardson-number-based diffusivites,
! as well as an additional ``entrainment'' diffusivity.
! The Louis diffusivities for momentum, $K_m$, and for heat
! and moisture, $K_h$, are defined at the interior layer edges. For LM layers,
! we define diffusivities at the base of the top LM-1 layers. All indexing
! is from top to bottom of the atmosphere.
!
!
! The Richardson number, Ri, is defined at the same edges as the diffusivities.
! $$
! {\rm Ri}_l = \frac{ \frac{g}{\left(\overline{\theta_v}\right)_l}\left(\frac{\delta \theta_v}{\delta z}\right)_l }
! { \left(\frac{\delta {\bf |V|}}{\delta z}\right)^2_l }, \, \, l=1,LM-1
! $$
! where $\theta_v=\theta(1+\epsilon q)$ is the virtual potential temperature,
! $\epsilon=\frac{M_a}{M_w}-1$, $M_a$ and $M_w$ are the molecular weights of
! dry air and water, and $q$ is the specific humidity.
! $\delta \theta_v$ is the difference of $\theta_v$ in the layers above and below the edge
! at which Ri$_l$ is defined; $\overline{\theta_v}$ is their average.
!
! The diffusivities at the layer edges have the form:
! $$
! K^m_l = (\ell^2_m)_l \left(\frac{\delta {\bf |V|}}{\delta z}\right)_l f_m({\rm Ri}_l)
! $$
! and
! $$
! K^h_l = (\ell^2_h)_l \left(\frac{\delta {\bf |V|}}{\delta z}\right)_l f_h({\rm Ri}_l),
! $$
! where $k$ is the Von Karman constant, and $\ell$ is the
! Blackdar(1962) length scale, also defined at the layer edges.
!
! Different turbulent length scales can be used for heat and momentum.
! in both cases, we use the traditional formulation:
! $$
! (\ell_{(m,h)})_l = \frac{kz_l}{1 + \frac{kz_l}{\lambda_{(m,h)}}},
! $$
! where, near the surface, the scale is proportional to $z_l$, the height above
! the surface of edge level $l$, and far from the surface it approaches $\lambda$.
! The length scale $\lambda$ is usually taken to be a constant (order 150 m), assuming
! the same scale for the outre boundary layer and the free atmosphere. We make it
! a function of height, reducing its value in the free atmosphere. The momentum
! length scale written as:
! $$
! \lambda_m = \max(\lambda_1 e^{\left(\frac{z_l}{z_T}\right)^2}, \lambda_2)
! $$
! where $\lambda_2 \le \lambda_1$ and $z_T$ is the top of the boundary layer.
! The length scale for heat and other scalers is taken as: $\lambda_h = \sqrt\frac{3d}{2} \lambda_m$,
! following the scheme used at ECMWF.
!
! The two universal functions of the Richardson number, $f_m$ and $f_h$,
! are taken from Louis et al (1982). For unstable conditions (Ri$\le 0$),
! they are:
! $$
! f_m = (1 - 2b \psi)
! $$
! and
! $$
! f_h = (1 - 3b \psi),
! $$
! where
! $$
! \psi = \frac{ {\rm Ri} }{ 1 + 3bC(z)\sqrt{-{\rm Ri}} },
! $$
! and
! $$
! C(z)=
! $$
! For stable condition (Ri$\ge 0$), they are
! $$
! f_m = \frac{1}{1.0 + \frac{2b{\rm Ri}}{\psi}}
! $$
! and
! $$
! f_h = \frac{1}{1.0 + 3b{\rm Ri}\psi},
! $$
! where
! $$
! \psi = \sqrt{1+d{\rm Ri}}.
! $$
! As in Louis et al (1982), the parameters appearing in these are taken
! as $b = c = d = 5$.
!EOP
! ErrLog Variables
character(len=ESMF_MAXSTR) :: IAm = "TurbulenceGetks"
integer :: STATUS
! Locals
real, dimension(size(KH,1),size(KH,2),size(KH,3)) :: ALH, ALM, DZ, DT, TM, PS, LAMBDAM_X, LAMBDAH_X
integer :: LM,L
real :: Zchoke
! Begin...
!===> Number of layers; edge levels will be one less (LM-1).
LM = size(ZZ,3)
!===> Quantities needed for Richardson number
DZ = (ZZ (:,:,1:LM-1) - ZZ (:,:,2:LM))
TM = (PV (:,:,1:LM-1) + PV (:,:,2:LM))*0.5
DT = (PV (:,:,1:LM-1) - PV (:,:,2:LM))
DU = (UU (:,:,1:LM-1) - UU (:,:,2:LM))**2 + &
(VV (:,:,1:LM-1) - VV (:,:,2:LM))**2
!===> Limits on distance between layer centers and vertical shear at edges.
DZ = max( DZ, MINTHICK)
DU = sqrt(DU)/DZ
!===> Richardson number ( RI = G*(DTheta_v/DZ) / (Theta_v*|DV/DZ|^2) )
RI = MAPL_GRAV*(DT/DZ)/(TM*( max(DU, MINSHEAR)**2))
!===> Blackadar(1962) length scale: $1/l = 1/(kz) + 1/\lambda$
LAMBDAM_X = MAX( LAMBDAM * EXP( -(ZE / ZKMENV )**2 ) , LAMBDAM2 )
LAMBDAH_X = MAX( LAMBDAH * EXP( -(ZE / ZKHENV )**2 ) , LAMBDAH2 )
ALM = ( MAPL_KARMAN*ZE/( 1.0 + MAPL_KARMAN*(ZE/LAMBDAM_X) ) )**2
ALH = ( MAPL_KARMAN*ZE/( 1.0 + MAPL_KARMAN*(ZE/LAMBDAH_X) ) )**2
if (associated( ALH_DIAG ) ) then
ALH_DIAG(:,:,0) = 0.0
ALH_DIAG(:,:,1:LM-1) = SQRT( ALH )
ALH_DIAG(:,:,LM) = 0.0
endif
!===> Unstable case: Uses (3.14, 3.18, 3.27) in Louis-scheme
! should approach (3.13) for small -RI.
where( RI < 0.0 )
PS = ( (ZZ(:,:,1:LM-1)/ZZ(:,:,2:LM))**(1./3.) - 1.0 ) ** 3
PS = ALH*sqrt( PS/(ZE*(DZ**3)) )
PS = RI/(1.0 + (3.0*LOUIS*LOUIS)*PS*sqrt(-RI))
KH = 1.0 - (LOUIS*3.0)*PS
KM = 1.0 - (LOUIS*2.0)*PS
endwhere
!===> Choke off unstable KH below Zchoke (m). JTB 2/2/06
Zchoke = 500.
where( (RI < 0.0) .and. (ZE < Zchoke ) )
KH = KH * (( ZE / Zchoke )**3)
endwhere
!===> Stable case
where( RI >= 0.0 )
PS = sqrt (1.0 + LOUIS *RI )
KH = 1.0 / (1.0 + (LOUIS*3.0)*RI*PS)
KM = PS / (PS + (LOUIS*2.0)*RI )
endwhere
!===> DIMENSIONALIZE Kz and LIMIT DIFFUSIVITY
ALM = DU*ALM
ALH = DU*ALH
KM = min(KM*ALM, KHMMAX)
KH = min(KH*ALH, KHMMAX)
RETURN_(ESMF_SUCCESS)
end subroutine LOUIS_KS
end module GEOS_TurbulenceGridCompMod