! +-======-+ ! Copyright (c) 2003-2007 United States Government as represented by ! the Admistrator of the National Aeronautics and Space Administration. ! All Rights Reserved. ! ! THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, ! REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN ! COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS ! REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). ! THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN ! INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR ! REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, ! DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED ! HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE ! RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. ! ! Government Agency: National Aeronautics and Space Administration ! Government Agency Original Software Designation: GSC-15354-1 ! Government Agency Original Software Title: GEOS-5 GCM Modeling Software ! User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov ! Government Agency Point of Contact for Original Software: ! Dale Hithon, SRA Assistant, (301) 286-2691 ! ! +-======-+ ! $Id$ #include "MAPL_Generic.h" !#define PRINT_STATES #define FULLPHYSICS #define GCM #define debug 0 ! Held-Suarez is not in module GEOSGCM????? #if defined HS #define GEOS_physicsGridCompMod GEOS_hsGridCompMod #undef FULLPHYSICS #endif #if defined SCM #define GEOS_superdynGridCompMod GEOS_singcolGridCompMod #undef GCM #endif !============================================================================= !BOP ! !MODULE: GEOS_AgcmGridCompMod -- A Module to combine Supedynamics and Physics Gridded Components ! !INTERFACE: module GEOS_AgcmGridCompMod ! !USES: use ESMF use MAPL_Mod use GEOS_TopoGetMod use GEOS_superdynGridCompMod, only: SDYN_SetServices => SetServices use GEOS_physicsGridCompMod, only: PHYS_SetServices => SetServices use MAPL_OrbGridCompMod, only: ORB_SetServices => SetServices use m_chars, only: uppercase use GEOS_RemapMod, only: myremap => remap implicit none private ! !PUBLIC MEMBER FUNCTIONS: public SetServices !============================================================================= ! !DESCRIPTION: This gridded component (GC) combines the Superdynamics GC, ! and Physics GC into a new composite Agcm GC. !\begin{verbatim} ! DUDT .... Mass-Weighted U-Wind Tendency (Pa m /s) ! DVDT .... Mass-Weighted V-Wind Tendency (Pa m /s) ! DPEDT ... Edge-Pressure Tendency (Pa /s) ! DTDT .... Mass-Weighted Temperature Tendency (Pa K /s) ! TRACER .. Friendly Tracers (unknown) !\end{verbatim} !EOP integer :: SDYN integer :: PHYS integer :: ORB type CONNECT_IAUcoeffs real, pointer :: dfi(:) => NULL() integer :: istep end type CONNECT_IAUcoeffs ! Wrapper for extracting internal state ! ------------------------------------- type IAU_coeffs type (CONNECT_IAUcoeffs), pointer :: PTR end type IAU_coeffs type CONNECT_ANAnBKG private type (MAPL_HorzTransform) :: ANA2BKG type (MAPL_HorzTransform) :: BKG2ANA type (ESMF_Grid) :: GRIDana type (ESMF_Alarm) :: AnaTendStartAlarm character(len=ESMF_MAXSTR) :: gridAnaName real, pointer :: phis_bkg(:,:) integer :: IM integer :: JM integer :: LM logical :: do_transforms=.false. logical :: initialized=.false. end type CONNECT_ANAnBKG ! Wrapper for extracting internal state ! ------------------------------------- type ANAnBKG type (CONNECT_ANAnBKG), pointer :: PTR end type ANAnBKG contains !BOP ! !IROUTINE: SetServices -- Sets ESMF services for this component ! !INTERFACE: subroutine SetServices ( GC, RC ) ! !ARGUMENTS: type(ESMF_GridComp), intent(INOUT) :: GC ! gridded component integer, optional :: RC ! return code ! !DESCRIPTION: The SetServices for the Physics GC needs to register its ! Initialize and Run. It uses the MAPL\_Generic construct for defining ! state specs and couplings among its children. In addition, it creates the ! children GCs (SURF, CHEM, RADIATION, MOIST, TURBULENCE) and runs their ! respective SetServices. !EOP !============================================================================= ! ! ErrLog Variables character(len=ESMF_MAXSTR) :: IAm integer :: STATUS character(len=ESMF_MAXSTR) :: COMP_NAME ! Locals integer :: I integer :: RST logical :: ANA_TS type (ESMF_Config) :: CF type (MAPL_MetaComp), pointer :: MAPL character(len=ESMF_MAXSTR) :: ReplayMode type (CONNECT_ANAnBKG), pointer :: upd_internal_state type (ANAnBKG) :: wrap type (CONNECT_IAUcoeffs), pointer :: iau_coeffs_internal_state type (IAU_coeffs) :: wrap_iau_coeffs !============================================================================= ! Begin... ! Get my name and set-up traceback handle ! --------------------------------------- call ESMF_GridCompGet( GC, NAME=COMP_NAME, RC=STATUS ) VERIFY_(STATUS) Iam = trim(COMP_NAME) // 'SetServices' ! Register services for this component ! ------------------------------------ call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_INITIALIZE, Initialize, RC=STATUS ) VERIFY_(STATUS) call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, Run , RC=STATUS ) VERIFY_(STATUS) ! Get the configuration from the component !----------------------------------------- call ESMF_GridCompGet( GC, CONFIG = CF, RC=STATUS ) VERIFY_(STATUS) ! Set the state variable specs. ! ----------------------------- call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, I, Label="ANALYZE_TS:", default=0, RC=STATUS) VERIFY_(STATUS) ANA_TS = I /= 0 if (ANA_TS) then RST = MAPL_RestartRequired else RST = MAPL_RestartSkip end if call MAPL_GetResource(MAPL, ReplayMode, Label='REPLAY_MODE:', default="NoReplay", RC=STATUS ) VERIFY_(STATUS) if(ANA_TS .and. ( adjustl(ReplayMode) /= "Exact" .and. & adjustl(ReplayMode) /= "Regular" ) ) then ASSERT_( adjustl(ReplayMode) == "NoReplay" ) endif !BOS ! !IMPORT STATE: call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'DUDT', & LONG_NAME = 'eastward_wind_analysis_increment', & UNITS = 'm s-1', & DIMS = MAPL_DimsHorzVert, & FIELD_TYPE = MAPL_VectorField, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'DVDT', & LONG_NAME = 'northward_wind_analysis_increment', & UNITS = 'm s-1', & DIMS = MAPL_DimsHorzVert, & FIELD_TYPE = MAPL_VectorField, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'DTDT', & LONG_NAME = 'temperature_analysis_increment', & UNITS = 'K', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'DPEDT', & LONG_NAME = 'edge_pressure_analysis_increment', & UNITS = 'Pa', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationEdge, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'DQVDT', & LONG_NAME = 'specific_humidity_analysis_increment', & UNITS = 'kg kg-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'DO3DT', & LONG_NAME = 'ozone_analysis_increment', & UNITS = 'mol mol-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'DTSDT', & LONG_NAME = 'skin_temperature_increment', & UNITS = 'K', & RESTART = RST, & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) ! !INTERNAL STATE: call MAPL_AddInternalSpec ( gc, & SHORT_NAME = 'DUDT', & LONG_NAME = 'eastward_wind_bias_tendency', & UNITS = 'm s-2', & FRIENDLYTO = trim(COMP_NAME), & DIMS = MAPL_DimsHorzVert, & FIELD_TYPE = MAPL_VectorField, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddInternalSpec ( gc, & SHORT_NAME = 'DVDT', & LONG_NAME = 'northward_wind_bias_tendency', & UNITS = 'm s-2', & FRIENDLYTO = trim(COMP_NAME), & DIMS = MAPL_DimsHorzVert, & FIELD_TYPE = MAPL_VectorField, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddInternalSpec ( gc, & SHORT_NAME = 'DTDT', & LONG_NAME = 'temperature_bias_tendency', & UNITS = 'K s-1', & FRIENDLYTO = trim(COMP_NAME), & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddInternalSpec ( gc, & SHORT_NAME = 'DPEDT', & LONG_NAME = 'edge_pressure_bias_tendency', & UNITS = 'Pa s-1', & FRIENDLYTO = trim(COMP_NAME), & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationEdge, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddInternalSpec ( gc, & SHORT_NAME = 'DQVDT', & LONG_NAME = 'specific_humidity_bias_tendency', & UNITS = 'kg kg-1 s-1', & FRIENDLYTO = trim(COMP_NAME), & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddInternalSpec ( gc, & SHORT_NAME = 'DO3DT', & LONG_NAME = 'ozone_bias_tendency', & UNITS = 'mol mol-1 s-1', & FRIENDLYTO = trim(COMP_NAME), & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddInternalSpec ( gc, & SHORT_NAME = 'DTSDT', & LONG_NAME = 'skin_temperature_tendency', & UNITS = 'K s-1', & FRIENDLYTO = trim(COMP_NAME), & RESTART = RST, & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) ! !EXPORT STATE: call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DPSDT_CONSTRAINT', & LONG_NAME = 'surface_pressure_adjustment_due_to_constraint', & UNITS = 'Pa s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DUDT_ANA', & LONG_NAME = 'total_eastward_wind_analysis_tendency', & UNITS = 'm s-2', & DIMS = MAPL_DimsHorzVert, & FIELD_TYPE = MAPL_VectorField, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DVDT_ANA', & LONG_NAME = 'total_northward_wind_analysis_tendency', & UNITS = 'm s-2', & DIMS = MAPL_DimsHorzVert, & FIELD_TYPE = MAPL_VectorField, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DTDT_ANA', & LONG_NAME = 'total_temperature_analysis_tendency', & UNITS = 'K s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DPEDT_ANA', & LONG_NAME = 'total_edge_pressure_analysis_tendency', & UNITS = 'Pa s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationEdge, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DQVDT_ANA', & LONG_NAME = 'total_specific_humidity_vapor_analysis_tendency', & UNITS = 'kg kg-1 s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DQLDT_ANA', & LONG_NAME = 'total_specific_humidity_liquid_analysis_tendency', & UNITS = 'kg kg-1 s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DQIDT_ANA', & LONG_NAME = 'total_specific_humidity_ice_analysis_tendency', & UNITS = 'kg kg-1 s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DO3DT_ANA', & LONG_NAME = 'total_ozone_analysis_tendency', & UNITS = 'mol mol-1 s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DTSDT_ANA', & LONG_NAME = 'total_skin_temperature_tendency', & UNITS = 'K s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DTHVDTFILINT', & LONG_NAME = 'vertically_integrated_thv_adjustment_from_filling',& UNITS = 'K kg m-2 s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'PERES', & LONG_NAME = 'vertically_integrated_cpt_tendency_residual', & UNITS = 'W m-2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'PEFILL', & LONG_NAME = 'vertically_integrated_cpt_adjustment_from_filling',& UNITS = 'W m-2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'QTFILL', & LONG_NAME = 'vertically_integrated_total_water_adjustment_from_filling', & UNITS = 'kg m-2 s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'QVFILL', & LONG_NAME = 'vertically_integrated_qv_adjustment_from_filling', & UNITS = 'kg m-2 s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'QLFILL', & LONG_NAME = 'vertically_integrated_ql_adjustment_from_filling', & UNITS = 'kg m-2 s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'QIFILL', & LONG_NAME = 'vertically_integrated_qi_adjustment_from_filling', & UNITS = 'kg m-2 s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'OXFILL', & LONG_NAME = 'vertically_integrated_ox_adjustment_from_filling', & UNITS = 'kg m-2 s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TROPP_EPV', & LONG_NAME = 'tropopause_pressure_based_on_EPV_estimate', & UNITS = 'Pa', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TROPP_THERMAL', & LONG_NAME = 'tropopause_pressure_based_on_thermal_estimate', & UNITS = 'Pa', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TROPP_BLENDED', & LONG_NAME = 'tropopause_pressure_based_on_blended_estimate', & UNITS = 'Pa', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TROPT', & LONG_NAME = 'tropopause_temperature_using_blended_TROPP_estimate', & UNITS = 'K', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TROPQ', & LONG_NAME = 'tropopause_specific_humidity_using_blended_TROPP_estimate', & UNITS = 'kg kg-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TQV', & LONG_NAME = 'total_precipitable_water_vapor', & UNITS = 'kg m-2' , & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TQI', & LONG_NAME = 'total_precipitable_ice_water', & UNITS = 'kg m-2' , & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TQL', & LONG_NAME = 'total_precipitable_liquid_water', & UNITS = 'kg m-2' , & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'TOX', & LONG_NAME = 'total_column_odd_oxygen', & UNITS = 'kg m-2' , & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'MASS', & LONG_NAME = 'atmospheric_mass', & UNITS = 'kg m-2' , & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'KE', & LONG_NAME = 'vertically_integrated_kinetic_energy', & UNITS = 'J m-2' , & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'CPT', & LONG_NAME = 'vertically_integrated_enthalpy', & UNITS = 'J m-2' , & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'THV', & LONG_NAME = 'vertically_integrated_virtual_potential_temperature', & UNITS = 'K' , & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'QLTOT', & LONG_NAME = 'mass_fraction_of_cloud_liquid_water', & UNITS = 'kg kg-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'QITOT', & LONG_NAME = 'mass_fraction_of_cloud_ice_water', & UNITS = 'kg kg-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'PHIS', & LONG_NAME = 'surface geopotential height', & UNITS = 'm+2 s-2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'SGH', & LONG_NAME = 'isotropic stdv of GWD topography', & UNITS = 'm', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'GWDVARX', & LONG_NAME = 'east-west variance of GWD topography', & UNITS = 'm+2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'GWDVARY', & LONG_NAME = 'north-south variance of GWD topography', & UNITS = 'm+2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'GWDVARXY', & LONG_NAME = 'SW-NE variance of GWD topography', & UNITS = 'm+2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'GWDVARYX', & LONG_NAME = 'NW-SE variance of GWD topography', & UNITS = 'm+2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'TRBVAR', & LONG_NAME = 'isotropic variance of TRB topography', & UNITS = 'm+2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'VARFLT', & LONG_NAME = 'isotropic variance of filtered topography', & UNITS = 'm+2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) ! Create childrens gridded components and invoke their SetServices ! ---------------------------------------------------------------- #ifdef SCM SDYN = MAPL_AddChild(GC, NAME='SCMDYNAMICS', SS=SDYN_SetServices, RC=STATUS) VERIFY_(STATUS) #else SDYN = MAPL_AddChild(GC, NAME='SUPERDYNAMICS', SS=SDYN_SetServices, RC=STATUS) VERIFY_(STATUS) #endif PHYS = MAPL_AddChild(GC, NAME='PHYSICS', SS=PHYS_SetServices, RC=STATUS) VERIFY_(STATUS) ORB = MAPL_AddChild(GC, NAME='ORBIT', SS=ORB_SetServices, RC=STATUS) VERIFY_(STATUS) ! Export for IAU or Analysis purposes ! ----------------------------------- ! call MAPL_AddExportSpec ( GC, & ! SHORT_NAME = 'PHIS', & ! CHILD_ID = SDYN, & ! RC=STATUS ) ! VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'AREA', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'AK', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'BK', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'PLE', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec( GC, & SHORT_NAME = 'PS', & CHILD_ID = SDYN, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'DELP', & CHILD_ID = SDYN, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'PE', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'PT', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'TV', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'T', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'U', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'V', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'U_DGRID', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'V_DGRID', & CHILD_ID = SDYN, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'O3PPMV', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'OX', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'Q', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'QCTOT', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'U10N', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'V10N', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'SNOMAS', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'WET1', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'TSOIL1', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'LWI', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'Z0', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'TS', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'TRANA', & CHILD_ID = PHYS, & RC = STATUS) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'FRLAND', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'FRLANDICE', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'FRLAKE', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'FROCEAN', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( GC, & SHORT_NAME = 'FRACI', & CHILD_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) !EOS ! Set internal connections between the childrens IMPORTS and EXPORTS ! ------------------------------------------------------------------ call MAPL_AddConnectivity ( GC, & SRC_NAME = (/'U ','V ','TH ','T ', & 'ZLE ','PS ','TA ','QA ', & 'SPEED ','DZ ','PLE ', & 'PREF ','TROPP_BLENDED','S ', & 'PV ','OMEGA '/), & DST_NAME = (/'U ','V ','TH ','T ', & 'ZLE ','PS ','TA ','QA ', & 'SPEED ','DZ ','PLE ', & 'PREF ','TROPP ','S ', & 'PV ','OMEGA '/), & DST_ID = PHYS, & SRC_ID = SDYN, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddConnectivity ( GC, & SRC_NAME = 'PLE', & DST_NAME = 'PLEINST', & SRC_ID = SDYN, & DST_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddConnectivity ( GC, SRC_NAME = 'AREA', DST_NAME = 'AREA', & SRC_ID = SDYN, DST_ID = PHYS, RC=STATUS ) VERIFY_(STATUS) ! Bundle of quantities to be advected !------------------------------------ call MAPL_AddConnectivity ( GC, & SRC_NAME = 'TRADV', & DST_NAME = 'TRADV', & SRC_ID = PHYS, & DST_ID = SDYN, & RC=STATUS ) VERIFY_(STATUS) ! Orbital Component Bundle call MAPL_AddConnectivity( GC, & SRC_NAME='SATORB', & DST_NAME='SATORB', & SRC_ID = ORB, & DST_ID = PHYS, & RC=STATUS ) VERIFY_(STATUS) #ifdef SCMSURF call MAPL_AddConnectivity ( GC, & SHORT_NAME = (/'TSKINOBS','QSKINOBS','LHOBS ','SHOBS '/), & DST_ID = PHYS, & SRC_ID = SDYN, & RC=STATUS ) VERIFY_(STATUS) #endif ! We Terminate these IMPORTS which are manually filled !----------------------------------------------------- call MAPL_TerminateImport ( GC, & SHORT_NAME = (/'DUDT ','DVDT ','DTDT ','DPEDT ','DQVANA','DQLANA','DQIANA','DOXANA','PHIS '/), & CHILD = SDYN, & RC=STATUS ) VERIFY_(STATUS) call MAPL_TerminateImport ( GC, & SHORT_NAME = (/'VARFLT','PHIS ','SGH ', 'DTSDT '/), & CHILD = PHYS, & RC=STATUS ) VERIFY_(STATUS) ! Allocate this instance of the internal state and put it in wrapper ! ------------------------------------------------------------------ allocate( upd_internal_state, stat=status ) VERIFY_(STATUS) wrap%ptr => upd_internal_state allocate( iau_coeffs_internal_state, stat=status ) VERIFY_(STATUS) wrap_iau_coeffs%ptr => iau_coeffs_internal_state ! Save pointer to the wrapped internal state in the GC ! ---------------------------------------------------- call ESMF_UserCompSetInternalState ( GC, 'UPD_STATE', wrap, status ) VERIFY_(STATUS) call ESMF_UserCompSetInternalState ( GC, 'IAU_COEFFS', wrap_iau_coeffs, status ) VERIFY_(STATUS) call MAPL_GenericSetServices ( GC, RC=STATUS ) VERIFY_(STATUS) ! Clocks !------- call MAPL_TimerAdd(GC, name="INITIALIZE" ,RC=STATUS) VERIFY_(STATUS) call MAPL_TimerAdd(GC, name="RUN" ,RC=STATUS) VERIFY_(STATUS) ! All done !--------- RETURN_(ESMF_SUCCESS) end subroutine SetServices !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !BOP ! !IROUTINE: Initialize -- Initialize method for the composite Agcm Gridded Component ! !INTERFACE: subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) ! !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: !EOP ! ErrLog Variables character(len=ESMF_MAXSTR) :: IAm integer :: STATUS character(len=ESMF_MAXSTR) :: COMP_NAME ! Local derived type aliases type (MAPL_MetaComp), pointer :: STATE type (ESMF_State), pointer :: GIM(:) type (ESMF_State), pointer :: GEX(:) type (ESMF_Field) :: FIELD type (ESMF_Time) :: CurrTime, RingTime type (ESMF_TimeInterval) :: TIMEINT type (ESMF_Alarm) :: ALARM type (ESMF_Alarm) :: ALARM4D type (ESMF_Config) :: cf integer :: I, NQ real :: POFFSET, DT real, pointer, dimension(:,:) :: PHIS,SGH,VARFLT,PTR real, pointer, dimension(:,:,:) :: TEND! character(len=ESMF_MAXSTR) :: replayMode real :: RPL_INTERVAL real :: RPL_SHUTOFF real :: IAU4dFREQ character(len=ESMF_MAXSTR), parameter :: INITIALIZED_EXPORTS(3) = & (/'PHIS ', 'SGH ', 'VARFLT' /) ! ============================================================================= ! Begin... ! Get the target components name and set-up traceback handle. ! ----------------------------------------------------------- call ESMF_GridCompGet ( GC, name=COMP_NAME, config=cf, RC=STATUS ) VERIFY_(STATUS) Iam = trim(COMP_NAME) // "Initialize" ! Get my MAPL_Generic state !-------------------------- call MAPL_GetObjectFromGC ( GC, STATE, RC=STATUS) VERIFY_(STATUS) call MAPL_TimerOn(STATE,"INITIALIZE") ! Call Initialize for every Child call MAPL_GenericInitialize ( GC, IMPORT, EXPORT, CLOCK, RC=STATUS) VERIFY_(STATUS) call MAPL_TimerOn(STATE,"TOTAL") ! Get children and their im/ex states from my generic state. !---------------------------------------------------------- call MAPL_Get ( STATE, GIM=GIM, GEX=GEX, RC=STATUS ) VERIFY_(STATUS) ! Make sure that the physics tendencies are allocated !---------------------------------------------------- call MAPL_GetPointer(GEX(PHYS), TEND, 'DUDT' , ALLOC=.true., rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(GEX(PHYS), TEND, 'DVDT' , ALLOC=.true., rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(GEX(PHYS), TEND, 'DTDT' , ALLOC=.true., rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(GEX(PHYS), TEND, 'DPEDT', ALLOC=.true., rc=STATUS) VERIFY_(STATUS) ! Fill Childrens TOPO variables and Diagnostics !---------------------------------------------- call MAPL_GetPointer(EXPORT, PHIS, 'PHIS', ALLOC=.true., rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, SGH, 'SGH', ALLOC=.true., rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, VARFLT, 'VARFLT', ALLOC=.true., rc=STATUS) VERIFY_(STATUS) ! PHIS ... !--------- call ESMF_StateGet( GIM(SDYN), 'PHIS', FIELD, rc=STATUS ) VERIFY_(STATUS) Call GEOS_TopoGet ( cf, MEAN=FIELD, rc=STATUS ) VERIFY_(STATUS) call ESMF_StateGet( GIM(PHYS), 'PHIS', FIELD, rc=STATUS ) VERIFY_(STATUS) Call GEOS_TopoGet ( cf, MEAN=FIELD, rc=STATUS ) VERIFY_(STATUS) call ESMF_FieldGet (FIELD, localDE=0, farrayPtr=PTR, rc = status) PHIS = PTR ! GWDVAR ... !----------- call ESMF_StateGet( GIM(PHYS), 'SGH', FIELD, rc=STATUS ) VERIFY_(STATUS) Call GEOS_TopoGet ( cf, GWDVAR=FIELD, rc=STATUS ) VERIFY_(STATUS) call ESMF_FieldGet (FIELD, localDE=0, farrayPtr=PTR, rc = status) SGH = PTR ! TRBVAR ... !----------- call ESMF_StateGet( GIM(PHYS), 'VARFLT', FIELD, rc=STATUS ) VERIFY_(STATUS) Call GEOS_TopoGet ( cf, TRBVAR=FIELD, rc=STATUS ) VERIFY_(STATUS) call ESMF_FieldGet (FIELD, localDE=0, farrayPtr=PTR, rc = status) VARFLT = PTR ! ====================================================================== !ALT: the next section addresses the problem when export variables have been ! assigned values during Initialize. To prevent "connected" exports ! being overwritten by DEFAULT in the Import spec in the other component ! we label them as being "initailized by restart". A better solution ! would be to move the computation to phase 2 of Initialize and ! eliminate this section alltogether ! ====================================================================== DO I = 1, size(INITIALIZED_EXPORTS) call ESMF_StateGet(EXPORT,INITIALIZED_EXPORTS(I), FIELD, RC=STATUS) VERIFY_(STATUS) call MAPL_AttributeSet(field, NAME="MAPL_InitStatus", & VALUE=MAPL_InitialRestart, RC=STATUS) VERIFY_(STATUS) END DO ! Initialize Predictor Alarm !--------------------------- call ESMF_ClockGet(clock, currTime=currTime, rc=status) VERIFY_(STATUS) call MAPL_GetResource( STATE, DT, Label="RUN_DT:", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource( STATE, POFFSET, Label="PREDICTOR_OFFSET:", default=21600. , RC=STATUS) VERIFY_(STATUS) call ESMF_TimeIntervalSet(TIMEINT, S=nint (POFFSET), RC=STATUS) VERIFY_(STATUS) ringTime = currTime+TIMEINT call ESMF_TimeIntervalSet(TIMEINT, S=nint(DT) , RC=STATUS) VERIFY_(STATUS) ALARM = ESMF_AlarmCreate( name='PredictorAlarm', & CLOCK = CLOCK, & RingInterval = TIMEINT , & RingTime = ringTime, & RC = STATUS ) VERIFY_(STATUS) if(ringTime == currTime) then call ESMF_AlarmRingerOn(Alarm, rc=status) VERIFY_(STATUS) end if call MAPL_StateAlarmAdd(STATE,ALARM,RC=status) VERIFY_(STATUS) call MAPL_GetResource(STATE, ReplayMode, 'REPLAY_MODE:', default="NoReplay", RC=STATUS ) VERIFY_(STATUS) if(adjustl(ReplayMode)=="Exact") then call MAPL_GetResource(STATE, RPL_SHUTOFF, 'REPLAY_SHUTOFF:', default=4000*21600., RC=STATUS ) VERIFY_(STATUS) call ESMF_TimeIntervalSet(TIMEINT, S=nint(RPL_SHUTOFF), RC=STATUS) VERIFY_(STATUS) ALARM = ESMF_AlarmCreate(name='ReplayShutOff', clock=CLOCK, & ringInterval=TIMEINT, sticky=.false., & RC=STATUS ) VERIFY_(STATUS) call MAPL_GetResource(STATE, RPL_INTERVAL, 'REPLAY_INTERVAL:', default=21600., RC=STATUS ) VERIFY_(STATUS) call ESMF_TimeIntervalSet(TIMEINT, S=nint(RPL_INTERVAL), RC=STATUS) VERIFY_(STATUS) ALARM = ESMF_AlarmCreate(name='ExactReplay', clock=CLOCK, & RingTime = currTime, & ringInterval=TIMEINT, sticky=.false., & RC=STATUS ) VERIFY_(STATUS) call ESMF_AlarmRingerOn(ALARM, rc=status) VERIFY_(STATUS) ASSERT_(POFFSET == RPL_INTERVAL) end if ! Create 4dIAU alarm ! ------------------ call ESMF_ClockGet(clock, currTime=currTime, rc=status) VERIFY_(STATUS) call MAPL_GetResource( STATE, IAU4dFREQ, Label="4DIAU_FREQUENCY:", default=3600. , RC=STATUS) VERIFY_(STATUS) call ESMF_TimeIntervalSet(TIMEINT, S=nint (IAU4dFREQ), RC=STATUS) VERIFY_(STATUS) ALARM4D = ESMF_AlarmCreate( name='4DIAUalarm', & CLOCK = CLOCK, & RingInterval = TIMEINT , & RingTime = currTime, & STICKY = .FALSE., & RC = STATUS ) VERIFY_(STATUS) call ESMF_AlarmRingerOn(Alarm4D, rc=status) VERIFY_(STATUS) call MAPL_TimerOff(STATE,"TOTAL") call MAPL_TimerOff(STATE,"INITIALIZE") #ifdef PRINT_STATES call WRITE_PARALLEL ( trim(Iam)//": IMPORT State" ) if ( MAPL_am_I_root() ) call ESMF_StatePrint ( IMPORT, rc=STATUS ) call WRITE_PARALLEL ( trim(Iam)//": EXPORT State" ) if ( MAPL_am_I_root() ) call ESMF_StatePrint ( EXPORT, rc=STATUS ) #endif RETURN_(ESMF_SUCCESS) end subroutine Initialize !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !BOP ! !IROUTINE: Run -- Run method for the composite Agcm Gridded Component ! !INTERFACE: subroutine Run ( GC, IMPORT, EXPORT, CLOCK, RC ) ! !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: !EOP ! ErrLog Variables character(len=ESMF_MAXSTR) :: IAm integer :: STATUS character(len=ESMF_MAXSTR) :: COMP_NAME ! Local derived type aliases type (MAPL_MetaComp), pointer :: STATE type (ESMF_GridComp), pointer :: GCS(:) type (ESMF_State), pointer :: GIM(:) type (ESMF_State), pointer :: GEX(:) type (ESMF_State) :: INTERNAL type (ESMF_Alarm) :: ALARM type (ESMF_Alarm) :: ALARM4D type (ESMF_TimeInterval) :: TINT type (ESMF_FieldBundle) :: Bundle type (ESMF_FieldBundle) :: Advect_Bundle type (ESMF_Grid) :: grid real, pointer, dimension(:) :: PREF => null() real, pointer, dimension(:,:,:) :: U => null() real, pointer, dimension(:,:,:) :: V => null() real, pointer, dimension(:,:,:) :: T => null() real, pointer, dimension(:,: ) :: TS => null() real, pointer, dimension(:,:,:) :: Q => null() real, pointer, dimension(:,:,:) :: QLLS => null() real, pointer, dimension(:,:,:) :: QLCN => null() real, pointer, dimension(:,:,:) :: QILS => null() real, pointer, dimension(:,:,:) :: QICN => null() real, pointer, dimension(:,:,:) :: PLE => null() real, pointer, dimension(:,:,:) :: EPV => null() real, pointer, dimension(:,:,:) :: QLTOT => null() real, pointer, dimension(:,:,:) :: QITOT => null() real, pointer, dimension(:,:,:) :: DPEDT => null() real, pointer, dimension(:,:,:) :: DTDT => null() real, pointer, dimension(:,:,:) :: TENDAN => null() #if 0 ! Simple Nudging Based on Predictor Methodology ! --------------------------------------------- real, pointer, dimension(:,:,:) :: u_ana => null() real, pointer, dimension(:,:,:) :: v_ana => null() real, pointer, dimension(:,:,:) :: t_ana => null() real, pointer, dimension(:,:,:) :: p_ana => null() real, pointer, dimension(:,:,:) :: q_ana => null() real, pointer, dimension(:,:,:) :: o3_ana => null() #endif real, allocatable, dimension(:,:) :: QFILL real, allocatable, dimension(:,:) :: QINT real, allocatable, dimension(:,:) :: QDP_BKG_INT real, allocatable, dimension(:,:) :: QDP_ANA_INT real*8, allocatable, dimension(:,:) :: SUMKE real*8, allocatable, dimension(:,:) :: SUMCPT1, SUMCPT2 real*8, allocatable, dimension(:,:) :: SUMTHV1, SUMTHV2 real*8, allocatable, dimension(:,:,:) :: PKE real*8, allocatable, dimension(:,:,:) :: PKZ real, pointer, dimension(:,:) :: AREA => null() real, pointer, dimension(:,:) :: OXFILL => null() real, pointer, dimension(:,:) :: QTFILL => null() real, pointer, dimension(:,:) :: QVFILL => null() real, pointer, dimension(:,:) :: QIFILL => null() real, pointer, dimension(:,:) :: QLFILL => null() real, pointer, dimension(:,:) :: TQV => null() real, pointer, dimension(:,:) :: TQI => null() real, pointer, dimension(:,:) :: TQL => null() real, pointer, dimension(:,:) :: TOX => null() real, pointer, dimension(:,:) :: TROPP1 => null() real, pointer, dimension(:,:) :: TROPP2 => null() real, pointer, dimension(:,:) :: TROPP3 => null() real, pointer, dimension(:,:) :: TROPT => null() real, pointer, dimension(:,:) :: TROPQ => null() real, pointer, dimension(:,:) :: MASS => null() real, pointer, dimension(:,:) :: KE => null() real, pointer, dimension(:,:) :: CPT => null() real, pointer, dimension(:,:) :: THV => null() real, pointer, dimension(:,:) :: PERES => null() real, pointer, dimension(:,:) :: PEFILL => null() real, pointer, dimension(:,:) :: PEPHY_SDYN => null() ! D(CpT)DT from ADD_INCS (SuperDYNamics) real, pointer, dimension(:,:) :: PEPHY_PHYS => null() ! D(CpT)DT from PHYSics real, pointer, dimension(:,:) :: DTHVDTFILINT => null() real, pointer, dimension(:,:) :: DTHVDTPHYINT => null() real, pointer, dimension(:,:) :: DQVDTPHYINT => null() real, pointer, dimension(:,:) :: DQLDTPHYINT => null() real, pointer, dimension(:,:) :: DQIDTPHYINT => null() real, pointer, dimension(:,:) :: DOXDTPHYINT => null() real, pointer, dimension(:,:,:) :: DP real, pointer, dimension(:,:,:) :: PL real, pointer, dimension(:,:,:) :: FC real, pointer, dimension(:,:,:) :: TROP real, pointer, dimension(:,:,:) :: XXINC real, pointer, dimension(:,:,:) :: DQVANA real, pointer, dimension(:,:,:) :: DQLANA real, pointer, dimension(:,:,:) :: DQIANA real, pointer, dimension(:,:,:) :: DOXANA real, pointer, dimension(:,:,:) :: DQIMPORT => null() real, pointer, dimension(:,:) :: ptr2d real, pointer, dimension(:,:,:) :: ptr3d real, pointer, dimension(:) :: AK real, pointer, dimension(:) :: BK real*8, allocatable, dimension(:,:) :: sumq real*8, allocatable, dimension(:,:) :: sum_qdp_bkg real*8, allocatable, dimension(:,:) :: sum_qdp_ana real, allocatable, dimension(:,:,:) :: qdp_bkg real, allocatable, dimension(:,:,:) :: qdp_ana real, allocatable, dimension(:,:,:) :: ple_ana real, allocatable, dimension(:,:,:) :: dp_ana real, allocatable, dimension(:,:,:) :: tdpold real, allocatable, dimension(:,:,:) :: tdpnew real*8 :: gamma real*8 :: qdp_bkg_ave real*8 :: qdp_ana_ave real*8 :: qint_ana_ave real*8 :: qint_bkg_ave real :: DT integer :: IM, JM, LM, L, TYPE, ISFCST integer :: NumFriendly integer :: K integer :: I logical :: DasMode logical :: DO_PREDICTOR logical :: LAST_CORRECTOR integer :: CONSTRAIN_DAS real :: ALPHA, BETA, TAUANL, DTX, IAUcoeff real :: ALPHAQ, BETAQ real :: ALPHAO, BETAO real :: ALF, BET real :: EPS character(len=ESMF_MAXSTR) :: ANA_IS_WEIGHTED logical :: IS_WEIGHTED character(len=ESMF_MAXSTR), pointer :: Names(:) character(len=ESMF_MAXSTR) :: STRING integer, parameter :: FREERUN = 0 integer, parameter :: PREDICTOR = 1 integer, parameter :: CORRECTOR = 2 integer, parameter :: FORECAST = 3 integer :: unit logical :: is_ringing logical,save :: is_shutoff=.false. character(len=ESMF_MAXSTR) :: FILENAME character(len=ESMF_MAXSTR) :: FileTmpl character(len=ESMF_MAXSTR) :: replayFile character(len=ESMF_MAXSTR) :: replayMode character(len=ESMF_MAXSTR) :: rplMode type(ESMF_Time) :: currTime CHARACTER(LEN=ESMF_MAXSTR) :: fieldName logical :: DO_4DIAU !============================================================================= ! Begin... ! Get the target components name and set-up traceback handle. ! ----------------------------------------------------------- Iam = "Run" call ESMF_GridCompGet ( GC, name=COMP_NAME, RC=STATUS ) VERIFY_(STATUS) Iam = trim(COMP_NAME) // trim(Iam) ! Get my MAPL_Generic state !-------------------------- call MAPL_GetObjectFromGC ( GC, STATE, RC=STATUS) VERIFY_(STATUS) call MAPL_TimerOn (STATE,"TOTAL") call MAPL_TimerOn (STATE,"RUN" ) ! Get children and their im/ex states from my generic state. !---------------------------------------------------------- call MAPL_Get ( STATE, GCS=GCS, GIM=GIM, GEX=GEX, & INTERNAL_ESMF_STATE=INTERNAL, & IM=IM, JM=JM, LM=LM, & RC=STATUS ) VERIFY_(STATUS) call ESMF_GridCompGet(GC, grid=grid, rc=status) VERIFY_(STATUS) ! Get the specific 4dIAU alarm !----------------------------- call ESMF_ClockGetAlarm(clock, alarmname='4DIAUalarm', alarm=Alarm4D, rc=status) VERIFY_(STATUS) ! Get the specific IAU alarm !--------------------------- call MAPL_StateAlarmGet(STATE, ALARM, NAME='PredictorAlarm', RC=STATUS) VERIFY_(STATUS) ! Set the various time scales !---------------------------- call MAPL_GetResource( STATE, DT, Label="RUN_DT:", RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, ALPHA, Label="ALPHA:", default=0.0, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, BETA, Label="BETA:", default=1.0, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, ALPHAQ, Label="ALPHAQ:", default=ALPHA, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, BETAQ, Label="BETAQ:", default=BETA, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, ALPHAO, Label="ALPHAO:", default=ALPHA, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, BETAO, Label="BETAO:", default=BETA, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, TAUANL, Label="TAUANL:", default=21600., RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, ISFCST, Label="IS_FCST:", default=0, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, CONSTRAIN_DAS, Label="CONSTRAIN_DAS:", default=1, RC=STATUS); VERIFY_(STATUS) call MAPL_GetResource( STATE, ANA_IS_WEIGHTED, Label="ANA_IS_WEIGHTED:", default='NO', RC=STATUS) VERIFY_(STATUS) ANA_IS_WEIGHTED = uppercase(ANA_IS_WEIGHTED) IS_WEIGHTED = adjustl(ANA_IS_WEIGHTED)=="YES" .or. adjustl(ANA_IS_WEIGHTED)=="NO" ASSERT_( IS_WEIGHTED ) IS_WEIGHTED = adjustl(ANA_IS_WEIGHTED)=="YES" call MAPL_GetResource( STATE, STRING, LABEL="IMPORT_RESTART_FILE:", RC=STATUS) IF (STATUS == ESMF_SUCCESS) THEN DasMode = .true. ELSE DasMode = .false. END IF call MAPL_GetResource( STATE, ReplayMode, 'REPLAY_MODE:', default="NoReplay", RC=STATUS ) VERIFY_(STATUS) ! NoReplay, Exact, Intermittent, Regular rplMode = adjustl(ReplayMode) if(rplMode=="Regular" .or. rplMode == "Exact") then DasMode = .true. end if ! Set type of update !------------------- if (ISFCST/=0 ) then TYPE = FORECAST else if(.not. DasMode) then TYPE = FREERUN else DO_PREDICTOR = ESMF_AlarmIsRinging( ALARM, rc=status) VERIFY_(STATUS) LAST_CORRECTOR = ESMF_AlarmWillRingNext( ALARM, rc=status) VERIFY_(STATUS) REPLAYING: if (rplMode == "Regular") then call ESMF_ClockGetAlarm(clock, 'startReplay', alarm, rc=status) VERIFY_(STATUS) LAST_CORRECTOR = ESMF_AlarmWillRingNext( ALARM, rc=status) VERIFY_(STATUS) else if (rplMode == "Exact") then call ESMF_AlarmRingerOff(ALARM, RC=STATUS) VERIFY_(STATUS) DO_PREDICTOR = .FALSE. call MAPL_GetResource ( STATE, FileTmpl,'REPLAY_FILE:', RC=STATUS ) VERIFY_(STATUS) ! If replay alarm is ringing, we need to reset state !--------------------------------------------------- if (is_shutoff) then ! once this alarm rings, is_shutoff will remain true for the rest of the run call MAPL_GetPointer(IMPORT,ptr3d,'DUDT',RC=STATUS) ptr3d=0.0 call MAPL_GetPointer(IMPORT,ptr3d,'DVDT',RC=STATUS) ptr3d=0.0 call MAPL_GetPointer(IMPORT,ptr3d,'DTDT',RC=STATUS) ptr3d=0.0 call MAPL_GetPointer(IMPORT,ptr3d,'DPEDT',RC=STATUS) ptr3d=0.0 call MAPL_GetPointer(IMPORT,ptr3d,'DQVDT',RC=STATUS) ptr3d=0.0 call MAPL_GetPointer(IMPORT,ptr3d,'DO3DT',RC=STATUS) ptr3d=0.0 call MAPL_GetPointer(IMPORT,ptr2d,'DTSDT',RC=STATUS) if(associated(ptr2d))then ptr2d=0.0 endif else call ESMF_ClockGetAlarm(Clock,'ReplayShutOff',Alarm,rc=Status) VERIFY_(status) is_shutoff = ESMF_AlarmWillRingNext( Alarm,rc=status ) VERIFY_(status) endif call ESMF_ClockGetAlarm(Clock,'ExactReplay',Alarm,rc=Status) VERIFY_(status) LAST_CORRECTOR = ESMF_AlarmWillRingNext( ALARM, rc=status) VERIFY_(STATUS) is_ringing = ESMF_AlarmIsRinging( Alarm,rc=status ) VERIFY_(status) is_ringing = is_ringing .and. (.not. is_shutoff) TIME_TO_REPLAY: if(is_ringing) then call ESMF_ClockGet(Clock, CurrTime=currTime, rc=Status) VERIFY_(status) ! when testing check alarms, predictor/corrector/last corrector ! template to file call MAPL_GetCurrentFile(FILETMPL=filetmpl, TIME=currTime, FILENAME=ReplayFile, & RC=STATUS) VERIFY_(status) unit = getfile(ReplayFile, FORM="unformatted", rc=status) VERIFY_(STATUS) ! read import from file call MAPL_VarRead(UNIT=UNIT, STATE=IMPORT, RC=STATUS) VERIFY_(STATUS) call FREE_FILE(unit, rc=status) VERIFY_(STATUS) end if TIME_TO_REPLAY end if REPLAYING if(DO_PREDICTOR) then TYPE = PREDICTOR else TYPE = CORRECTOR end if end if ! Get Names Associated with Friendly Analysis Bundle !--------------------------------------------------- call ESMF_StateGet(GEX(PHYS), 'TRANA', BUNDLE, RC=STATUS ) VERIFY_(STATUS) call ESMF_FieldBundleGet(BUNDLE,FieldCount=NumFriendly, RC=STATUS) VERIFY_(STATUS) ASSERT_(NumFriendly==2) allocate(Names(NumFriendly), stat=STATUS) VERIFY_(STATUS) call ESMF_FieldBundleGet(BUNDLE, fieldNameList=Names, RC=STATUS) VERIFY_(STATUS) ! Prepare for update !------------------- call MAPL_GetPointer( GEX(SDYN), PREF,'PREF',rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), PLE, 'PLE', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), U, 'U', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), V, 'V', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), T, 'T', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), AK, 'AK', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), BK, 'BK', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(PHYS), TS, 'TS', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), PEPHY_SDYN, 'PEPHY', alloc=.true., rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(PHYS), PEPHY_PHYS, 'PEPHY', alloc=.true., rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), DTHVDTPHYINT, 'DTHVDTPHYINT', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(PHYS), DQVDTPHYINT, 'DQVDTPHYINT', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(PHYS), DQLDTPHYINT, 'DQLDTPHYINT', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(PHYS), DQIDTPHYINT, 'DQIDTPHYINT', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(PHYS), DOXDTPHYINT, 'DOXDTPHYINT', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), AREA, 'AREA', rc=STATUS ) VERIFY_(STATUS) allocate( PL(IM,JM,LM),STAT=STATUS ) VERIFY_(STATUS) allocate( DP(IM,JM,LM),STAT=STATUS ) VERIFY_(STATUS) PL = 0.5*(PLE(:,:,1:LM)+PLE(:,:,0:LM-1)) DP = PLE(:,:,1:LM)-PLE(:,:,0:LM-1) ALF = ALPHA BET = BETA ! If so, overwrite increments from analysis by recreating them on the fly ! ----------------------------------------------------------------------- call get_iau_coeff(IAUcoeff) DO_4DIAU = ESMF_AlarmIsRinging( ALARM4D, rc=status) VERIFY_(STATUS) if(DO_4DIAU .and. (TYPE == CORRECTOR) ) then call ESMF_AlarmRingerOff(ALARM4D, RC=STATUS) VERIFY_(STATUS) call update_ainc_(RC=STATUS) VERIFY_(STATUS) endif #if 0 ! Simple Nudging Based on Predictor Methodology ! --------------------------------------------- if(TYPE == CORRECTOR ) then call MAPL_GetPointer( GEX(SDYN), PLE, 'PLE', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), U, 'U', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), V, 'V', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), T, 'T', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer(IMPORT,ptr3d,'DUDT',RC=STATUS) VERIFY_(STATUS) allocate( u_ana(size(ptr3d,1),size(ptr3d,2),size(ptr3d,3)) ) u_ana = ptr3d ptr3d = ptr3d - U call MAPL_GetPointer(IMPORT,ptr3d,'DVDT',RC=STATUS) VERIFY_(STATUS) allocate( v_ana(size(ptr3d,1),size(ptr3d,2),size(ptr3d,3)) ) v_ana = ptr3d ptr3d = ptr3d - V call MAPL_GetPointer(IMPORT,ptr3d,'DTDT',RC=STATUS) VERIFY_(STATUS) allocate( t_ana(size(ptr3d,1),size(ptr3d,2),size(ptr3d,3)) ) t_ana = ptr3d ptr3d = ptr3d - T call MAPL_GetPointer(IMPORT,ptr3d,'DPEDT',RC=STATUS) VERIFY_(STATUS) allocate( p_ana(size(ptr3d,1),size(ptr3d,2),size(ptr3d,3)) ) p_ana = ptr3d ptr3d = ptr3d - PLE do K=1,NumFriendly NULLIFY(Q) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) STRING = TRIM(Names(K)) fieldName = MAPL_RmQualifier(STRING) if(TRIM(fieldName) == 'OX') then call MAPL_GetPointer(IMPORT,ptr3d,'DO3DT',RC=STATUS) VERIFY_(STATUS) allocate( o3_ana(size(ptr3d,1),size(ptr3d,2),size(ptr3d,3)) ) o3_ana = ptr3d ptr3d = ptr3d - Q*1.e6 endif if(NAMES(K)=='Q') then call MAPL_GetPointer(IMPORT,ptr3d,'DQVDT',RC=STATUS) VERIFY_(STATUS) allocate( q_ana(size(ptr3d,1),size(ptr3d,2),size(ptr3d,3)) ) q_ana = ptr3d ptr3d = ptr3d - Q endif enddo endif #endif ! Load Analysis Increments into Imports for Physics and Dynamics State Variables !------------------------------------------------------------------------------- call DO_UPDATE_ANA2D ('DTSDT', PHYS) ! Note: Mass-Weighting is done for DTDT, No Mass-Weighting for DUDT,DVDT,DPEDT ! ---------------------------------------------------------------------------- call DO_UPDATE_ANA3D ('DUDT' , SDYN, PREF) call DO_UPDATE_ANA3D ('DVDT' , SDYN, PREF) call DO_UPDATE_ANA3D ('DPEDT', SDYN, PREF, CONSTRAIN_DAS = CONSTRAIN_DAS) call MAPL_GetPointer(GIM(SDYN), DPEDT, 'DPEDT', rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(GIM(SDYN), DTDT, 'DTDT', rc=STATUS) VERIFY_(STATUS) call DO_UPDATE_ANA3D ('DTDT' , SDYN, PREF) if( is_weighted ) then allocate( tdpold( IM,JM, LM),STAT=STATUS ) ; VERIFY_(STATUS) allocate( tdpnew( IM,JM, LM),STAT=STATUS ) ; VERIFY_(STATUS) allocate( dp_ana( IM,JM, LM),STAT=STATUS ) ; VERIFY_(STATUS) allocate( ple_ana( IM,JM,0:LM),STAT=STATUS ) ; VERIFY_(STATUS) ! Create Proxies for Updated Pressure and Temperature due to Analysis !-------------------------------------------------------------------- ple_ana = ple + dt*dpedt dp_ana = ple_ana(:,:,1:LM)-ple_ana(:,:,0:LM-1) tdpold = T*DP tdpnew = ( T + dt*dtdt )*dp_ana dtdt = ( tdpnew - tdpold )/dt deallocate( tdpold ) deallocate( tdpnew ) deallocate( dp_ana ) deallocate( ple_ana ) endif ! Add Analysis Increment Directly to Friendlies !---------------------------------------------- if(TYPE /= FREERUN) then allocate(FC (IM,JM,LM),STAT=STATUS) VERIFY_(STATUS) do K=1,NumFriendly FC = 1.0 NULLIFY(Q) !ALT: ESMF requires that the data pointer is not associated call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) STRING = TRIM(Names(K)) fieldName = MAPL_RmQualifier(STRING) if(TRIM(fieldName) == 'OX') then ! PCHEM OX or GOCART::OX, for example. ALF = ALPHAO BET = BETAO DTX = DT*1.0E-6 ! Uncomment damping if problem with ozone ! --------------------------------------- ! do L=1,LM ! where(PL(:,:,L) < 100.0 .and. PL(:,:,L) > 0.0 ) ! FC(:,:,L) = exp(-1.5*(log10(PL(:,:,L))-2.0)**2) ! end where ! end do call MAPL_GetPointer(GIM(SDYN), DOXANA, 'DOXANA', rc=STATUS) VERIFY_(STATUS) DOXANA = Q call DO_Friendly (Q,'DO3DT',PREF) DOXANA = Q - DOXANA else ALF = ALPHAQ BET = BETAQ DTX = DT if(NAMES(K)=='Q') then ! Q ! Initialize DQVANA Diagnostics with Background QV !------------------------------------------------- call MAPL_GetPointer(GIM(SDYN), DQVANA, 'DQVANA', rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(GIM(SDYN), DQLANA, 'DQLANA', rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(GIM(SDYN), DQIANA, 'DQIANA', rc=STATUS) VERIFY_(STATUS) ! Get Pointers to QL & QI from Friendly Advection Bundle !------------------------------------------------------- call ESMF_StateGet(GEX(PHYS), 'TRADV', Advect_BUNDLE, RC=STATUS ) VERIFY_(STATUS) call ESMFL_BundleGetPointerToData( Advect_BUNDLE, 'QLLS', QLLS, RC=STATUS ) VERIFY_(STATUS) call ESMFL_BundleGetPointerToData( Advect_BUNDLE, 'QLCN', QLCN, RC=STATUS ) VERIFY_(STATUS) call ESMFL_BundleGetPointerToData( Advect_BUNDLE, 'QILS', QILS, RC=STATUS ) VERIFY_(STATUS) call ESMFL_BundleGetPointerToData( Advect_BUNDLE, 'QICN', QICN, RC=STATUS ) VERIFY_(STATUS) DQVANA = Q DQLANA = QLLS + QLCN DQIANA = QILS + QICN if(TYPE == CORRECTOR) then IF( CONSTRAIN_DAS == 1 ) then ! --------------------------- ! Create BKG Water Mass ! --------------------- allocate( qdp_bkg( IM,JM,LM ),STAT=STATUS ) ; VERIFY_(STATUS) do L=1,lm qdp_bkg(:,:,L) = ( q(:,:,L)+qlls(:,:,L)+qlcn(:,:,L)+qils(:,:,L)+qicn(:,:,L) )*dp(:,:,L) enddo ENDIF #if debug allocate( qint( IM,JM ),STAT=STATUS ) ; VERIFY_(STATUS) allocate( sumq( IM,JM ),STAT=STATUS ) ; VERIFY_(STATUS) sumq = 0.0_8 do L=1,lm sumq = sumq + ( q(:,:,L)+qlls(:,:,L)+qlcn(:,:,L)+qils(:,:,L)+qicn(:,:,L) )*dp(:,:,L) enddo qint = sumq call MAPL_AreaMean( qint_bkg_ave, qint, area, grid, rc=STATUS ) VERIFY_(STATUS) #endif ENDIF ! End CORRECTOR Test call DO_Friendly (Q,'DQVDT',PREF) if(TYPE == CORRECTOR) then IF( CONSTRAIN_DAS == 1 ) then ! --------------------------- ! Create Proxies for Updated Pressure and Temperature due to Analysis !-------------------------------------------------------------------- allocate( dp_ana(IM,JM, LM),STAT=STATUS ) ; VERIFY_(STATUS) allocate( ple_ana(IM,JM,0:LM),STAT=STATUS ) ; VERIFY_(STATUS) ple_ana = ple + dt*dpedt dp_ana = ple_ana(:,:,1:LM)-ple_ana(:,:,0:LM-1) ! Create ANA Water Mass ! --------------------- allocate( qdp_ana(IM,JM,LM),STAT=STATUS ); VERIFY_(STATUS) do L=1,lm qdp_ana(:,:,L) = ( q(:,:,L)+qlls(:,:,L)+qlcn(:,:,L)+qils(:,:,L)+qicn(:,:,L) )*dp_ana(:,:,L) enddo ! Vertically Integrate ANA & BKG Water Mass where they Differ ! ----------------------------------------------------------- allocate( sum_qdp_bkg( IM,JM ),STAT=STATUS ) ; VERIFY_(STATUS) allocate( sum_qdp_ana( IM,JM ),STAT=STATUS ) ; VERIFY_(STATUS) sum_qdp_bkg = 0.0_8 sum_qdp_ana = 0.0_8 do L=1,lm where( qdp_ana(:,:,L).ne.qdp_bkg(:,:,L) ) sum_qdp_bkg = sum_qdp_bkg + qdp_bkg(:,:,L) sum_qdp_ana = sum_qdp_ana + qdp_ana(:,:,L) end where enddo ! Compute Area-Mean Vertically Integrated BKG Water Mass ! ------------------------------------------------------ allocate( qdp_bkg_int( IM,JM ),STAT=STATUS ) ; VERIFY_(STATUS) where( sum_qdp_bkg.ne.0.0_8 ) qdp_bkg_int = sum_qdp_bkg elsewhere qdp_bkg_int = MAPL_UNDEF end where call MAPL_AreaMean( qdp_bkg_ave, qdp_bkg_int, area, grid, rc=STATUS ) VERIFY_(STATUS) ! Compute Area-Mean Vertically Integrated ANA Water Mass ! ------------------------------------------------------ allocate( qdp_ana_int( IM,JM ),STAT=STATUS ) ; VERIFY_(STATUS) where( sum_qdp_ana.ne.0.0_8 ) qdp_ana_int = sum_qdp_ana elsewhere qdp_ana_int = MAPL_UNDEF end where call MAPL_AreaMean( qdp_ana_ave, qdp_ana_int, area, grid, rc=STATUS ) VERIFY_(STATUS) ! Compute Dry-Mass Scaling Parameter ! ---------------------------------- if( real(qdp_bkg_ave,kind=4).ne.MAPL_UNDEF .and. & real(qdp_ana_ave,kind=4).ne.MAPL_UNDEF ) then gamma = real( qdp_bkg_ave / qdp_ana_ave , kind=4 ) else gamma = 1.0_8 endif ! Scale ANA Water Variables for Dry-Mass Conservation ! --------------------------------------------------- do L=1,lm where( qdp_ana(:,:,L).ne.qdp_bkg(:,:,L) ) q (:,:,L) = q (:,:,L) * gamma qlls(:,:,L) = qlls(:,:,L) * gamma qlcn(:,:,L) = qlcn(:,:,L) * gamma qils(:,:,L) = qils(:,:,L) * gamma qicn(:,:,L) = qicn(:,:,L) * gamma end where enddo #if debug sumq = 0.0_8 do L=1,lm sumq = sumq + ( q(:,:,L)+qlls(:,:,L)+qlcn(:,:,L)+qils(:,:,L)+qicn(:,:,L) )*dp_ana(:,:,L) enddo qint = sumq call MAPL_AreaMean( qint_ana_ave, qint, area, grid, rc=STATUS ) VERIFY_(STATUS) if(MAPL_AM_I_ROOT() ) then write(6,1001) qint_ana_ave,qint_bkg_ave,qint_ana_ave-qint_bkg_ave,gamma 1001 format(5x,'Global_QDP_ANA: ',g,' Global_QDP_BKG: ',g,' DIFF: ',g,' GAMMA: ',g) endif deallocate( qint ) deallocate( sumq ) #endif deallocate( dp_ana ) deallocate( ple_ana ) deallocate( qdp_bkg ) deallocate( qdp_ana ) deallocate( qdp_bkg_int ) deallocate( qdp_ana_int ) deallocate( sum_qdp_bkg ) deallocate( sum_qdp_ana ) ENDIF ! End CONSTRAIN_DAS Test ENDIF ! End CORRECTOR Test DQVANA = Q - DQVANA DQLANA = QLLS + QLCN - DQLANA DQIANA = QILS + QICN - DQIANA ! Update Tendency Diagnostic due to CONSTRAINTS ! --------------------------------------------- call MAPL_GetPointer ( EXPORT, TENDAN, 'DQVDT_ANA', rc=STATUS ) VERIFY_(STATUS) if(associated(TENDAN)) TENDAN = DQVANA/DT call MAPL_GetPointer ( EXPORT, TENDAN, 'DQLDT_ANA', rc=STATUS ) VERIFY_(STATUS) if(associated(TENDAN)) TENDAN = DQLANA/DT call MAPL_GetPointer ( EXPORT, TENDAN, 'DQIDT_ANA', rc=STATUS ) VERIFY_(STATUS) if(associated(TENDAN)) TENDAN = DQIANA/DT else call DO_Friendly (Q,'D'//trim(Names(K))//'DT',PREF) end if end if end do deallocate(FC) end if ! not free-running ! Make Sure EPV is Allocated for TROPOPAUSE Diagnostics !------------------------------------------------------ call MAPL_GetPointer ( EXPORT, TROPP1, 'TROPP_THERMAL', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, TROPP2, 'TROPP_EPV' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, TROPP3, 'TROPP_BLENDED', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, TROPT, 'TROPT', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, TROPQ, 'TROPQ', rc=STATUS ) VERIFY_(STATUS) if( associated(TROPP1) .or. & associated(TROPP2) .or. & associated(TROPP3) .or. & associated(TROPT) .or. & associated(TROPQ) ) then call MAPL_GetPointer( GEX(SDYN),EPV,'EPV',ALLOC=.true.,rc=STATUS ) VERIFY_(STATUS) endif ! Call basic run phase for both Child !------------------------------------- ! Call run for satellite orbits !------------------------------ call MAPL_TimerOn (STATE,"ORBIT" ) call ESMF_GridCompRun(GCS(ORB), importState=GIM(ORB), exportState=GEX(ORB), clock=CLOCK, userRC=STATUS) VERIFY_(STATUS) call MAPL_TimerOff(STATE,"ORBIT" ) call MAPL_TimerOn (STATE,"SUPERDYNAMICS" ) call ESMF_GridCompRun(GCS(SDYN), importState=GIM(SDYN), exportState=GEX(SDYN), clock=CLOCK, PHASE=1, userRC=STATUS) VERIFY_(STATUS) call MAPL_TimerOFF (STATE,"SUPERDYNAMICS" ) call MAPL_TimerOn (STATE,"PHYSICS" ) call ESMF_GridCompRun(GCS(PHYS), importState=GIM(PHYS), exportState=GEX(PHYS), clock=CLOCK, PHASE=1, userRC=STATUS) VERIFY_(STATUS) call MAPL_TimerOff(STATE,"PHYSICS" ) ! Load Physics Tendencies into Imports for RUN2 of Dynamics (ADD_INCS) !--------------------------------------------------------------------- call DO_UPDATE_PHY ('DUDT' ) call DO_UPDATE_PHY ('DVDT' ) call DO_UPDATE_PHY ('DPEDT') call DO_UPDATE_PHY ('DTDT' ) ! Run RUN2 of SuperDynamics (ADD_INCS) to add Physics Diabatic Tendencies !------------------------------------------------------------------------ call MAPL_TimerOn (STATE,"SUPERDYNAMICS" ) call ESMF_GridCompRun(GCS(SDYN), importState=GIM(SDYN), exportState=GEX(SDYN), clock=CLOCK, PHASE=2, userRC=STATUS) VERIFY_(STATUS) call MAPL_TimerOff(STATE,"SUPERDYNAMICS" ) ! Get Names Associated with Friendly Advection Bundle for Final Check for Negative Tracers !----------------------------------------------------------------------------------------- call ESMF_StateGet(GEX(PHYS), 'TRADV', BUNDLE, RC=STATUS ) VERIFY_(STATUS) call ESMF_FieldBundleGet(BUNDLE,FieldCount=NumFriendly, RC=STATUS) VERIFY_(STATUS) deallocate(Names) allocate(Names(NumFriendly), stat=STATUS) VERIFY_(STATUS) call ESMF_FieldBundleGet(BUNDLE, fieldNameList=Names, RC=STATUS) VERIFY_(STATUS) ! Get Pointers to Exports !------------------------ call MAPL_GetPointer ( EXPORT, QTFILL, 'QTFILL', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, QVFILL, 'QVFILL', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, QIFILL, 'QIFILL', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, QLFILL, 'QLFILL', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, OXFILL, 'OXFILL', rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, TOX , 'TOX' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, TQV , 'TQV' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, TQI , 'TQI' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, TQL , 'TQL' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, QLTOT , 'QLTOT' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, QITOT , 'QITOT' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, PERES , 'PERES' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, PEFILL , 'PEFILL' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, DTHVDTFILINT, 'DTHVDTFILINT' , rc=STATUS ) VERIFY_(STATUS) if(associated(QTFILL)) QTFILL = 0.0 if(associated(QIFILL)) QIFILL = 0.0 if(associated(QLFILL)) QLFILL = 0.0 if(associated(TQI) ) TQI = 0.0 if(associated(TQL) ) TQL = 0.0 if(associated(QLTOT) ) QLTOT = 0.0 if(associated(QITOT) ) QITOT = 0.0 allocate(QFILL(IM,JM) ,STAT=STATUS ) VERIFY_(STATUS) allocate(QINT (IM,JM) ,STAT=STATUS ) VERIFY_(STATUS) allocate( PKE(IM,JM,0:LM),STAT=STATUS ) VERIFY_(STATUS) allocate( PKZ(IM,JM,1:LM),STAT=STATUS ) VERIFY_(STATUS) PL = 0.5*(PLE(:,:,1:LM)+PLE(:,:,0:LM-1)) ! Recompute Updated Pressure DP = PLE(:,:,1:LM)-PLE(:,:,0:LM-1) ! Recompute Updated Pressure Thickness PKE = PLE**MAPL_KAPPA do L=1,LM PKZ(:,:,L) = ( PKE(:,:,L)-PKE(:,:,L-1) ) / ( MAPL_KAPPA*( log(PLE(:,:,L))-log(PLE(:,:,L-1)) ) ) enddo ! Initialize Vertically Integrated Values of CPT and THV (before QFILL updates) ! ----------------------------------------------------------------------------- if( associated(PEPHY_SDYN) .or. associated(PEFILL) ) then EPS = MAPL_RVAP/MAPL_RGAS-1.0 do K=1,NumFriendly NULLIFY(Q) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) if(NAMES(K)=='Q') then allocate( SUMCPT1(IM,JM),STAT=STATUS ) VERIFY_(STATUS) SUMCPT1 = 0.0 do L=1,LM SUMCPT1 = SUMCPT1 + MAPL_CP*T(:,:,L)*(1.0+EPS*Q(:,:,L))*DP(:,:,L) enddo exit endif enddo endif if( associated(DTHVDTPHYINT) .or. associated(DTHVDTFILINT) ) then EPS = MAPL_RVAP/MAPL_RGAS-1.0 do K=1,NumFriendly NULLIFY(Q) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) if(NAMES(K)=='Q') then allocate( SUMTHV1(IM,JM),STAT=STATUS ) VERIFY_(STATUS) SUMTHV1 = 0.0 do L=1,LM SUMTHV1 = SUMTHV1 + T(:,:,L)/PKZ(:,:,L)*(1.0+EPS*Q(:,:,L))*DP(:,:,L) enddo exit endif enddo endif ! Perform Final Check for Negative Friendlies ! ------------------------------------------- do K=1,NumFriendly NULLIFY(Q) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) STRING = TRIM(Names(K)) fieldName = MAPL_RmQualifier(STRING) ! Water Vapor ! ----------- if(NAMES(K)=='Q') then call FILL_Friendly ( Q,DP,QFILL,QINT ) if(associated(QVFILL)) QVFILL = QFILL if(associated(QTFILL)) QTFILL = QTFILL + QFILL if(associated(DQVDTPHYINT)) DQVDTPHYINT = DQVDTPHYINT + QFILL if(associated(TQV)) TQV = QINT endif ! Ice Water ! --------- if(NAMES(K)=='QICN') then call FILL_Friendly ( Q,DP,QFILL,QINT ) if(associated(QIFILL)) QIFILL = QIFILL + QFILL if(associated(QTFILL)) QTFILL = QTFILL + QFILL if(associated(DQIDTPHYINT)) DQIDTPHYINT = DQIDTPHYINT + QFILL if(associated(TQI)) TQI = TQI + QINT if(associated(QITOT)) QITOT = QITOT + Q endif if(NAMES(K)=='QILS') then call FILL_Friendly ( Q,DP,QFILL,QINT ) if(associated(QIFILL)) QIFILL = QIFILL + QFILL if(associated(QTFILL)) QTFILL = QTFILL + QFILL if(associated(DQIDTPHYINT)) DQIDTPHYINT = DQIDTPHYINT + QFILL if(associated(TQI)) TQI = TQI + QINT if(associated(QITOT)) QITOT = QITOT + Q endif ! Liquid Water ! ------------ if(NAMES(K)=='QLCN') then call FILL_Friendly ( Q,DP,QFILL,QINT ) if(associated(QLFILL)) QLFILL = QLFILL + QFILL if(associated(QTFILL)) QTFILL = QTFILL + QFILL if(associated(DQLDTPHYINT)) DQLDTPHYINT = DQLDTPHYINT + QFILL if(associated(TQL)) TQL = TQL + QINT if(associated(QLTOT)) QLTOT = QLTOT + Q endif if(NAMES(K)=='QLLS') then call FILL_Friendly ( Q,DP,QFILL,QINT ) if(associated(QLFILL)) QLFILL = QLFILL + QFILL if(associated(QTFILL)) QTFILL = QTFILL + QFILL if(associated(DQLDTPHYINT)) DQLDTPHYINT = DQLDTPHYINT + QFILL if(associated(TQL)) TQL = TQL + QINT if(associated(QLTOT)) QLTOT = QLTOT + Q endif ! Total Odd-Oxygen ! ---------------- if(TRIM(fieldName) == 'OX') then call FILL_Friendly ( Q,DP,QFILL,QINT ) if(associated(OXFILL)) OXFILL = QFILL *(MAPL_O3MW/MAPL_AIRMW) if(associated(DOXDTPHYINT)) DOXDTPHYINT = DOXDTPHYINT + QFILL*(MAPL_O3MW/MAPL_AIRMW) if(associated(TOX)) TOX = QINT *(MAPL_O3MW/MAPL_AIRMW) endif enddo deallocate(QFILL) deallocate(QINT ) ! Compute Additional Diagnostics !------------------------------- call MAPL_GetPointer ( EXPORT, MASS , 'MASS' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, KE , 'KE' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, CPT , 'CPT' , rc=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer ( EXPORT, THV , 'THV' , rc=STATUS ) VERIFY_(STATUS) if( associated(MASS) ) MASS = (PLE(:,:,LM)-PLE(:,:,0)) * (1.0/MAPL_GRAV) if( associated(KE) ) then allocate( SUMKE(IM,JM),STAT=STATUS ) VERIFY_(STATUS) SUMKE = 0.0 do L=1,LM SUMKE = SUMKE + 0.5*( U(:,:,L)**2 + V(:,:,L)**2 )*DP(:,:,L) enddo KE = SUMKE * (1.0/MAPL_GRAV) deallocate(SUMKE) endif ! Instantaneous Values of CPT and THV are done here to include possible QFILL updates ! ----------------------------------------------------------------------------------- if( associated(CPT) .or. associated(PEPHY_SDYN) .or. associated(PEFILL) ) then EPS = MAPL_RVAP/MAPL_RGAS-1.0 do K=1,NumFriendly NULLIFY(Q) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) if(NAMES(K)=='Q') then allocate( SUMCPT2(IM,JM),STAT=STATUS ) VERIFY_(STATUS) SUMCPT2 = 0.0 do L=1,LM SUMCPT2 = SUMCPT2 + MAPL_CP*T(:,:,L)*(1.0+EPS*Q(:,:,L))*DP(:,:,L) enddo if( associated(CPT) ) CPT = SUMCPT2 * (1.0/MAPL_GRAV) if( associated(PEFILL) .or. & associated(PEPHY_SDYN) ) then SUMCPT1 = ( SUMCPT2-SUMCPT1 ) / (DT*MAPL_GRAV) if( associated(PEFILL) ) PEFILL = SUMCPT1 if( associated(PEPHY_SDYN) ) PEPHY_SDYN = PEPHY_SDYN + SUMCPT1 deallocate(SUMCPT1) endif deallocate(SUMCPT2) exit endif enddo endif if( associated(PERES) .and. & associated(PEPHY_SDYN) .and. & associated(PEPHY_PHYS) ) PERES = PEPHY_SDYN - PEPHY_PHYS if( associated(THV) .or. associated(DTHVDTPHYINT) .or. associated(DTHVDTFILINT) ) then EPS = MAPL_RVAP/MAPL_RGAS-1.0 do K=1,NumFriendly NULLIFY(Q) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) if(NAMES(K)=='Q') then allocate( SUMTHV2(IM,JM),STAT=STATUS ) VERIFY_(STATUS) SUMTHV2 = 0.0 do L=1,LM SUMTHV2 = SUMTHV2 + T(:,:,L)/PKZ(:,:,L)*(1.0+EPS*Q(:,:,L))*DP(:,:,L) enddo if( associated(THV) ) THV = SUMTHV2 * (MAPL_P00**MAPL_KAPPA) * (1.0/MAPL_GRAV) if( associated(DTHVDTFILINT) .or. & associated(DTHVDTPHYINT) ) then SUMTHV1 = ( SUMTHV2-SUMTHV1 ) * (MAPL_P00**MAPL_KAPPA) / (DT*MAPL_GRAV) if( associated(DTHVDTFILINT) ) DTHVDTFILINT = SUMTHV1 if( associated(DTHVDTPHYINT) ) DTHVDTPHYINT = DTHVDTPHYINT + SUMTHV1 deallocate(SUMTHV1) endif deallocate(SUMTHV2) exit endif enddo endif ! Compute Tropopause Diagnostics ! ------------------------------ if( associated(TROPP1) .or. & associated(TROPP2) .or. & associated(TROPP3) .or. & associated(TROPT) .or. & associated(TROPQ) ) then do K=1,NumFriendly NULLIFY(Q) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) if(NAMES(K)=='Q') then allocate( TROP(IM,JM,5),STAT=STATUS ) VERIFY_(STATUS) call tropovars ( IM,JM,LM,PLE,PL,T,Q,EPV,TROP(:,:,1),TROP(:,:,2),TROP(:,:,3),TROP(:,:,4),TROP(:,:,5) ) if( associated(TROPP1) ) TROPP1(:,:) = TROP(:,:,1) if( associated(TROPP2) ) TROPP2(:,:) = TROP(:,:,2) if( associated(TROPP3) ) TROPP3(:,:) = TROP(:,:,3) if( associated(TROPT ) ) TROPT (:,:) = TROP(:,:,4) if( associated(TROPQ ) ) TROPQ (:,:) = TROP(:,:,5) deallocate(TROP ) exit endif enddo endif ! Done !----- #if 0 ! Simple Nudging Based on Predictor Methodology ! --------------------------------------------- if(TYPE == CORRECTOR ) then call MAPL_GetPointer(IMPORT,ptr3d,'DUDT',RC=STATUS) VERIFY_(STATUS) ptr3d = u_ana call MAPL_GetPointer(IMPORT,ptr3d,'DVDT',RC=STATUS) VERIFY_(STATUS) ptr3d = v_ana call MAPL_GetPointer(IMPORT,ptr3d,'DTDT',RC=STATUS) VERIFY_(STATUS) ptr3d = t_ana call MAPL_GetPointer(IMPORT,ptr3d,'DPEDT',RC=STATUS) VERIFY_(STATUS) ptr3d = p_ana do K=1,NumFriendly NULLIFY(Q) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), Q, RC=STATUS ) VERIFY_(STATUS) STRING = TRIM(Names(K)) fieldName = MAPL_RmQualifier(STRING) if(TRIM(fieldName) == 'OX') then call MAPL_GetPointer(IMPORT,ptr3d,'DO3DT',RC=STATUS) VERIFY_(STATUS) ptr3d = o3_ana endif if(NAMES(K)=='Q') then call MAPL_GetPointer(IMPORT,ptr3d,'DQVDT',RC=STATUS) VERIFY_(STATUS) ptr3d = q_ana endif enddo deallocate( u_ana ) deallocate( v_ana ) deallocate( t_ana ) deallocate( p_ana ) deallocate( q_ana ) deallocate( o3_ana ) endif #endif deallocate(Names) deallocate(PL) deallocate(DP) deallocate(PKE) deallocate(PKZ) call MAPL_TimerOff(STATE,"RUN" ) call MAPL_TimerOff(STATE,"TOTAL") RETURN_(ESMF_SUCCESS) contains subroutine DO_UPDATE_ANA3D(NAME, COMP, PREF, CONSTRAIN_DAS) character*(*), intent(IN) :: NAME integer, intent(IN) :: COMP real, pointer, intent(IN) :: PREF(:) integer, optional, intent(IN) :: CONSTRAIN_DAS real, pointer, dimension(:,:) :: DPSDT_CONSTRAINT => null() real, pointer, dimension(:,:,:) :: TENDSD => null() real, pointer, dimension(:,:,:) :: TENDPH => null() real, pointer, dimension(:,:,:) :: TENDBS => null() real, pointer, dimension(:,:,:) :: TENDAN => null() real, pointer, dimension(:,:,:) :: ANAINC => null() real, allocatable, dimension(:,:,:) :: TENDANAL real, allocatable, dimension(:,:) :: dummy real, allocatable, dimension(:) :: ALFZ, BETZ real*8 :: qave1 real*8 :: qave2 real*8 :: qave3 integer :: L,LL,LU call MAPL_GetPointer(GIM(COMP), TENDSD, trim(NAME) , rc=STATUS) VERIFY_(STATUS) select case (TYPE) case (FREERUN) TENDSD = 0.0 case (PREDICTOR) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) TENDSD = TENDBS case (FORECAST) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) TENDBS = BET * TENDBS TENDSD = TENDBS case (CORRECTOR) LL = lbound(TENDSD,3) LU = ubound(TENDSD,3) allocate(TENDANAL(size(TENDSD,1),size(TENDSD,2),LL:LU), stat=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(IMPORT , ANAINC, trim(NAME), rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, DPSDT_CONSTRAINT, 'DPSDT_CONSTRAINT', rc=STATUS) VERIFY_(STATUS) TENDANAL = ANAINC*IAUcoeff ! No Constraints IF( present(CONSTRAIN_DAS) ) then IF(CONSTRAIN_DAS == 1 ) then ! -------------------------- allocate(dummy(size(TENDSD,1),size(TENDSD,2)), stat=STATUS) VERIFY_(STATUS) ! Method to Simply Re-Scale ! ------------------------- where( ANAINC(:,:,LU).ne.0.0 ) dummy = ANAINC(:,:,LU) ! ANAINC = PS_ANA - PS_BKG elsewhere dummy = MAPL_UNDEF endwhere call MAPL_AreaMean( qave1, dummy, area, grid, rc=STATUS ) ! qave1 = AreaMean( ANAINC ) VERIFY_(STATUS) where( ANAINC(:,:,LU).ne.0.0 ) dummy = PLE(:,:,LU) ! P_n elsewhere dummy = MAPL_UNDEF endwhere call MAPL_AreaMean( qave2, dummy, area, grid, rc=STATUS ) ! qave2 = AreaMean( P_n ) if( qave1.ne.MAPL_UNDEF .and. qave2.ne.MAPL_UNDEF ) then qave3 = qave2 + qave1*dt*IAUcoeff ! qave3 = AreaMean( P_n+1 = P_n + ANAINC*dt/tau ) else qave3 = MAPL_UNDEF endif if( qave3.ne.MAPL_UNDEF ) then where( ANAINC(:,:,LU).ne.0.0 ) dummy = ANAINC(:,:,LU)*real(qave2/qave3,kind=4) - dummy * real(qave1/qave3,kind=4) elsewhere dummy = 0.0 endwhere else dummy = 0.0 endif if(associated(DPSDT_CONSTRAINT)) then DPSDT_CONSTRAINT = ( dummy - ANAINC(:,:,LU) )*IAUcoeff endif DO L=LL,LU TENDANAL(:,:,L) = dummy(:,:)*BK(L) ENDDO TENDANAL = TENDANAL*IAUcoeff deallocate(dummy) ELSE if(associated(DPSDT_CONSTRAINT)) then DPSDT_CONSTRAINT = 0.0 endif ENDIF ENDIF TENDSD = (TENDBS + TENDANAL) if (LAST_CORRECTOR) then allocate( ALFZ(LL:LU), stat=STATUS) VERIFY_(STATUS) allocate( BETZ(LL:LU), stat=STATUS) VERIFY_(STATUS) DO L=LL,LU if( PREF(L).GT.10000.0 ) then ALFZ(L) = ALF ! PREF > 100 mb BETZ(L) = BET ! PREF > 100 mb elseif( PREF(L).LT.1000.0 ) then ALFZ(L) = 0.0 ! PREF < 10 mb BETZ(L) = 0.0 ! PREF < 10 mb else ALFZ(L) = ALF*( PREF(L)-1000.0 )/9000.0 BETZ(L) = BET*( PREF(L)-1000.0 )/9000.0 endif ENDDO DO L=LL,LU TENDBS(:,:,L) = BETZ(L)*TENDBS(:,:,L) + ALFZ(L)*TENDANAL(:,:,L) ENDDO deallocate(ALFZ) deallocate(BETZ) end if deallocate(TENDANAL) case default ASSERT_(.false.) end select ! Fill Total Increment Tendency (Current + Bias) Diagnostic ! --------------------------------------------------------- call MAPL_GetPointer ( EXPORT, TENDAN, trim(NAME)//'_ANA', rc=STATUS ) VERIFY_(STATUS) if(associated(TENDAN)) then TENDAN = TENDSD endif end subroutine DO_UPDATE_ANA3D subroutine DO_UPDATE_ANA2D(NAME, COMP) character*(*), intent(IN) :: NAME integer, intent(IN) :: COMP real, pointer, dimension(:,:) :: TENDSD => null() real, pointer, dimension(:,:) :: TENDPH => null() real, pointer, dimension(:,:) :: TENDBS => null() real, pointer, dimension(:,:) :: TENDAN => null() real, pointer, dimension(:,:) :: ANAINC => null() real, allocatable, dimension(:,:) :: TENDANAL call MAPL_GetPointer(GIM(COMP), TENDSD, trim(NAME) , rc=STATUS) VERIFY_(STATUS) select case (TYPE) case (FREERUN) TENDSD = 0.0 case (PREDICTOR) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) TENDSD = TENDBS case (FORECAST) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) TENDBS = BET * TENDBS TENDSD = TENDBS case (CORRECTOR) allocate( TENDANAL(size(TENDSD,1),size(TENDSD,2)), stat=STATUS ) VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(IMPORT , ANAINC, trim(NAME), rc=STATUS) VERIFY_(STATUS) TENDANAL = ANAINC*IAUcoeff TENDSD = (TENDBS + TENDANAL) if (LAST_CORRECTOR) then TENDBS = BET*TENDBS + ALF*TENDANAL end if deallocate(TENDANAL) case default ASSERT_(.false.) end select ! Fill Total Increment Tendency (Current + Bias) Diagnostic ! --------------------------------------------------------- call MAPL_GetPointer ( EXPORT, TENDAN, trim(NAME)//'_ANA', rc=STATUS ) VERIFY_(STATUS) if(associated(TENDAN)) then TENDAN = TENDSD endif end subroutine DO_UPDATE_ANA2D subroutine DO_UPDATE_PHY (NAME) character*(*), intent(IN) :: NAME real, pointer, dimension(:,:,:) :: TENDSD => null() real, pointer, dimension(:,:,:) :: TENDPH => null() call MAPL_GetPointer(GIM(SDYN), TENDSD, trim(NAME) , rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(GEX(PHYS), TENDPH, trim(NAME) , rc=STATUS) VERIFY_(STATUS) TENDSD = TENDPH end subroutine DO_UPDATE_PHY subroutine DO_Friendly(Q, NAME, PREF) character*(*), intent(IN ) :: NAME real, pointer, intent(IN) :: PREF(:) real, intent(INOUT) :: Q(:,:,:) real, pointer, dimension(:,:,:) :: TENDBS => null() real, pointer, dimension(:,:,:) :: TENDAN => null() real, pointer, dimension(:,:,:) :: ANAINC => null() real, allocatable, dimension(:,:,:) :: TENDANAL real, allocatable, dimension(:,:,:) :: QOLD real, allocatable, dimension(:) :: ALFZ, BETZ allocate( QOLD(IM,JM,LM), stat=STATUS) VERIFY_(STATUS) QOLD = Q ! Initialize Old Value for Total Tendency Diagnostic select case (TYPE) case (PREDICTOR) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) Q = Q + max( DTX*TENDBS*FC, -Q ) ! Prevent Negative Q case (FORECAST) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) TENDBS = BET * TENDBS Q = Q + max( DTX*TENDBS*FC, -Q ) ! Prevent Negative Q case (CORRECTOR) allocate(TENDANAL(IM,JM,LM), stat=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL , TENDBS, trim(NAME), rc=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(IMPORT , ANAINC, trim(NAME), rc=STATUS) VERIFY_(STATUS) TENDANAL = ANAINC*IAUcoeff Q = Q + max( DTX*(TENDBS + TENDANAL)*FC, -Q ) ! Prevent Negative Q if (LAST_CORRECTOR) then allocate( ALFZ(LM), stat=STATUS) VERIFY_(STATUS) allocate( BETZ(LM), stat=STATUS) VERIFY_(STATUS) DO L=1,LM if( PREF(L).GT.10000.0 ) then ALFZ(L) = ALF ! PREF > 100 mb BETZ(L) = BET ! PREF > 100 mb elseif( PREF(L).LT.1000.0 ) then ALFZ(L) = 0.0 ! PREF < 10 mb BETZ(L) = 0.0 ! PREF < 10 mb else ALFZ(L) = ALF*( PREF(L)-1000.0 )/9000.0 BETZ(L) = BET*( PREF(L)-1000.0 )/9000.0 endif ENDDO DO L=1,LM TENDBS(:,:,L) = BETZ(L)*TENDBS(:,:,L) + ALFZ(L)*TENDANAL(:,:,L) ENDDO deallocate(ALFZ) deallocate(BETZ) end if deallocate(TENDANAL) case default ASSERT_(.false.) end select ! Fill Total Increment Tendency (Current + Bias) Diagnostic ! --------------------------------------------------------- call MAPL_GetPointer ( EXPORT, TENDAN, trim(NAME)//'_ANA', rc=STATUS ) VERIFY_(STATUS) if(associated(TENDAN)) TENDAN = (Q-QOLD)/DT deallocate(QOLD) end subroutine DO_FRIENDLY subroutine FILL_Friendly ( Q,DP,QFILL,QINT ) real, intent(INOUT) :: Q(:,:,:) real, intent(IN ) :: DP(:,:,:) real, intent(OUT ) :: QFILL(:,:) real, intent(OUT ) :: QINT (:,:) real*8, allocatable, dimension(:,:) :: QTEMP1 real*8, allocatable, dimension(:,:) :: QTEMP2 allocate(QTEMP1(IM,JM), stat=STATUS) VERIFY_(STATUS) allocate(QTEMP2(IM,JM), stat=STATUS) VERIFY_(STATUS) QTEMP1 = 0.0 do L=1,LM QTEMP1(:,:) = QTEMP1(:,:) + Q(:,:,L)*DP(:,:,L) enddo where( Q < 0.0 ) Q = 0.0 QTEMP2 = 0.0 do L=1,LM QTEMP2(:,:) = QTEMP2(:,:) + Q(:,:,L)*DP(:,:,L) enddo where( qtemp2.ne.0.0_8 ) qtemp2 = max( qtemp1/qtemp2, 0.0_8 ) end where do L=1,LM Q(:,:,L) = Q(:,:,L)*qtemp2(:,:) enddo QTEMP2 = 0.0 do L=1,LM QTEMP2(:,:) = QTEMP2(:,:) + Q(:,:,L)*DP(:,:,L) enddo WHERE( QTEMP1 >= 0.0 ) QFILL = 0.0 ELSEWHERE QFILL = -QTEMP1 / (DT*MAPL_GRAV) END WHERE QINT = QTEMP2 / MAPL_GRAV deallocate(QTEMP1) deallocate(QTEMP2) end subroutine FILL_FRIENDLY subroutine get_iau_coeff(TNDCoeff) use m_chars, only: uppercase implicit none real,intent(OUT) :: TNDCoeff integer :: nsteps real :: POFFSET character(len=ESMF_MAXSTR) :: FilterType type (CONNECT_IAUcoeffs), pointer :: myCoeffs => NULL() type (IAU_coeffs) :: wrap call MAPL_GetResource( STATE, POFFSET, Label="PREDICTOR_OFFSET:", default=21600. , RC=STATUS) VERIFY_(STATUS) nsteps=nint(POFFSET/DT)+1 call MAPL_GetResource( STATE, FilterType, Label="FILTER_TYPE:", default="IAU" , RC=STATUS) VERIFY_(STATUS) FilterType = uppercase(FilterType) TNDcoeff = 1.0/TAUANL ! Get my internal private state. This contains the transforms ! between the background grid and the ANA grid, as well as ANA grid ! itself. !------------------------------------------------------------- call ESMF_UserCompGetInternalState(GC, 'IAU_COEFFS', wrap, status) VERIFY_(STATUS) myCoeffs => wrap%ptr ! Define coefficients to scale increments with ! NOTE: this is placed here to all 3dIAU to use this feature as well ! ------------------------------------------------------------------ if (.not.associated(myCoeffs%dfi)) then allocate(myCoeffs%dfi(nsteps)) myCoeffs%istep=0 call dfi_coeffs (FilterType,DT,POFFSET,TAUANL,nsteps,myCoeffs%dfi) if (MAPL_am_I_root()) then print*, 'DFI initialized for this many steps: ', nsteps do i=1,nsteps print*, 'i,LH-dfi-coeff=', i, myCoeffs%dfi(i),' (',1.0/(myCoeffs%dfi(i)*3600),' hrs)' enddo endif endif if (trim(FilterType)/="IAU") then myCoeffs%istep=myCoeffs%istep+1 if(myCoeffs%istep<=nsteps) then TNDcoeff=myCoeffs%dfi(myCoeffs%istep) else TNDcoeff=0.0 endif endif end subroutine get_iau_coeff subroutine update_ainc_(RC) use ESMF_CFIOMOD, only: ESMF_CFIOstrTemplate use ESMF_CFIOFileMod use GEOS_UtilsMod use GEOS_RemapMod, only: myremap => remap use m_chars, only: uppercase implicit none integer,optional, intent(OUT) :: RC ! Brief description: This routine converts either analysis increments or ! analysis fields to model increments corresponding to what IAU ! tendencies require. ! In case of handling analysis full fields: ! - convert background fields corresponding to analysis fields to ! analysis resolution. ! - calculate (IAU) increment on analysis grid ! - if needed, interpolate increments to model grid ! In case of handling analysis increment: ! - if needed, interpolate analysis increment to model grid ! - convert analysis increment to IAU increment (e.g., tv to td) ! The following name declarations will be easily unwired ... character(len=*), parameter :: incnames(7) = (/ 'sphu ', & 'u ', & 'v ', & 'tv ', & 'ozone', & 'ps ', & 'ts ' /) character(len=*), parameter :: ananames(9) = (/ 'sphu ', & 'u ', & 'v ', & 'tv ', & 'ozone', & 'delp ', & 'phis ', & 'ps ', & 'ts ' /) integer rank,ni,nt,nvars,natts integer IMana_World,JMana_World integer IMbkg_World,JMbkg_World integer IMana,JMana,LMana integer IMbkg,JMbkg,LMbkg integer NX,NY integer nymd,nhms integer vm_comm integer :: DIMS(ESMF_MAXGRIDDIM) integer :: K,NQ,FID integer :: idum integer :: method integer :: AINC_TIME(6) logical do_transforms logical l_cube logical l_nudge logical l_remap logical l_store_transforms logical fromANA2BKG logical IuseTS logical l_windfix type(ESMF_Config) :: CF type(ESMF_Field) :: Field type(ESMF_FieldBundle) :: RBUNDLE type(ESMF_Time) :: AincTime type(ESMF_VM) :: VM real,allocatable, dimension(:,:) :: dps real,allocatable, dimension(:,:,:) :: dp_inc real,allocatable, dimension(:,:,:) :: uptr real,allocatable, dimension(:,:,:) :: vptr ! real, pointer, dimension(:,:,:) :: qdum1 real, pointer, dimension(:,:,:) :: qdum2 ! real, pointer, dimension(:,:,:) :: u_bkg => NULL() real, pointer, dimension(:,:,:) :: u_ana => NULL() real, pointer, dimension(:,:,:) :: du_inc => NULL() real, pointer, dimension(:,:,:) :: v_bkg => NULL() real, pointer, dimension(:,:,:) :: v_ana => NULL() real, pointer, dimension(:,:,:) :: dv_inc => NULL() real, pointer, dimension(:,:,:) :: tv_bkg => NULL() real, pointer, dimension(:,:,:) :: tv_ana => NULL() real, pointer, dimension(:,:,:) :: dt_inc => NULL() real, pointer, dimension(:,:,:) :: q_bkg => NULL() real, pointer, dimension(:,:,:) :: q_ana => NULL() real, pointer, dimension(:,:,:) :: dq_inc => NULL() real, pointer, dimension(:,:,:) :: o3_bkg => NULL() real, pointer, dimension(:,:,:) :: o3_ana => NULL() real, pointer, dimension(:,:,:) :: do3_inc => NULL() real, pointer, dimension(:,:,:) :: ple_bkg => NULL() real, pointer, dimension(:,:,:) ::delp_bkg => NULL() real, pointer, dimension(:,:,:) :: ple_ana => NULL() real, pointer, dimension(:,:,:) ::dple_inc => NULL() real, pointer, dimension(:,:,:) :: pk_ana => NULL() real, pointer, dimension(:,:,:) :: pke_ana => NULL() real, pointer, dimension(:,:,:) :: pk_bkg => NULL() real, pointer, dimension(:,:,:) :: pke_bkg => NULL() real, pointer, dimension(:,:,:) :: dpk => NULL() real, pointer, dimension(:,:,:) :: thv_ana => NULL() real, pointer, dimension(:,:) :: ts_bkg => NULL() real, pointer, dimension(:,:) :: ts_ana => NULL() real, pointer, dimension(:,:) :: dts_inc => NULL() real, pointer, dimension(:,:) :: ps_bkg => NULL() real, pointer, dimension(:,:) ::phis_bkg real, pointer, dimension(:,:) ::phis_ana => NULL() ! real, pointer, dimension(:,:,:) :: dpkz => NULL() ! real, pointer, dimension(:,:) :: dts_aux => NULL() real, pointer, dimension(:,:,:) :: du_aux => NULL() real, pointer, dimension(:,:,:) :: dv_aux => NULL() real, pointer, dimension(:,:,:) :: dt_aux => NULL() real, pointer, dimension(:,:,:) :: dp_aux => NULL() real, pointer, dimension(:,:,:) :: dq_aux => NULL() real, pointer, dimension(:,:,:) :: do3_aux => NULL() real, pointer, dimension(:,:,:) :: dple_aux => NULL() ! real,pointer :: aptr2d(:,:) => NULL() ! analysis increment pointers real,pointer :: aptr3d(:,:,:)=> NULL() ! analysis increment pointers real,pointer :: gptr2d(:,:) => NULL() ! gcm background pointers real,pointer :: gptr3d(:,:,:)=> NULL() ! gcm background pointers ! real, allocatable, dimension(:,:) :: vintdiva real, allocatable, dimension(:,:) :: vintdivb real, allocatable, dimension(:,:) :: vintdivc ! character(len=ESMF_MAXSTR) :: name character(len=ESMF_MAXSTR) :: imstr, jmstr character(len=ESMF_MAXSTR) :: FILETMPL,AINCFILE character(len=ESMF_MAXSTR) :: STRING character(len=ESMF_MAXSTR) :: DATE character(len=ESMF_MAXSTR) :: NUDGE character(len=ESMF_MAXSTR) :: NUDGE_REMAP character(len=ESMF_MAXSTR) :: NUDGE_WINDFIX character(len=ESMF_MAXSTR) :: NUDGE_STORE_TRANSFORMS character(len=ESMF_MAXSTR) :: SKIP_FIRST type(ESMF_TimeInterval) :: Frequency type (MAPL_HorzTransform),pointer :: ANA2BKG => NULL() type (MAPL_HorzTransform),pointer :: BKG2ANA => NULL() type (CONNECT_ANAnBKG), pointer :: myANA => NULL() type (ANAnBKG) :: wrap type(ESMF_Grid), EXTERNAL :: AppGridCreateF integer :: ifirst data ifirst /0/ ! When assimilation period is over, do not even bother ... ! -------------------------------------------------------- call MAPL_GetResource(STATE, FILETMPL, LABEL="AINC_FILE:", default="NULL", RC=STATUS) VERIFY_(STATUS) ! If analysis (or increment) file not specified, move on ... ! ---------------------------------------------------------- if( trim(FILETMPL)=="NULL") return ! Get pointers to analysis tendencies ! ----------------------------------- call MAPL_GetPointer(import, du_inc, 'DUDT', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, dv_inc, 'DVDT', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, dt_inc, 'DTDT', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, dq_inc, 'DQVDT', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, do3_inc, 'DO3DT', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, dple_inc, 'DPEDT', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, dts_inc, 'DTSDT', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(STATE, SKIP_FIRST, LABEL="SKIP_FIRST:", default="NO", RC=STATUS) VERIFY_(STATUS) if(trim(SKIP_FIRST)=="YES") then ! when trying to "replay" to trajectory of free model ! either have agcm_import filled w/ zeros ! or bypass zero tendencies out in first ! pass if(ifirst<1)then du_inc =0.0 dv_inc =0.0 dt_inc =0.0 dq_inc =0.0 do3_inc =0.0 dple_inc=0.0 dts_inc =0.0 ifirst =ifirst+1 return endif else ! typicall, will expect tendency at initial time ! to be meaningful so, use what is in agcm_import ! in the first pass, from then on analysis should ! be read in at desired frequency if(ifirst<1)then ifirst =ifirst+1 return endif endif ! Get my internal private state. This contains the transforms ! between the background grid and the ANA grid, as well as ANA grid ! itself. !------------------------------------------------------------- call ESMF_UserCompGetInternalState(GC, 'UPD_STATE', wrap, status) VERIFY_(STATUS) myANA => wrap%ptr call MAPL_GetResource(STATE, NUDGE, LABEL="NUDGE_STATE:", default="NO", RC=STATUS) VERIFY_(STATUS) NUDGE = uppercase(NUDGE) l_nudge=trim(NUDGE)=="YES" call MAPL_GetResource(STATE, NUDGE_REMAP, LABEL="NUDGE_REMAP:", default="YES", RC=STATUS) VERIFY_(STATUS) NUDGE_REMAP = uppercase(NUDGE_REMAP) l_remap=trim(NUDGE_REMAP)=="YES" call MAPL_GetResource(STATE, NUDGE_WINDFIX, LABEL="NUDGE_WINDFIX:", default="YES", RC=STATUS) VERIFY_(STATUS) NUDGE_WINDFIX = uppercase(NUDGE_WINDFIX) l_windfix=trim(NUDGE_WINDFIX)=="YES" call MAPL_GetResource(STATE, NUDGE_STORE_TRANSFORMS, LABEL="NUDGE_STORE_TRANSFORMS:", default="YES", RC=STATUS) VERIFY_(STATUS) NUDGE_STORE_TRANSFORMS = uppercase(NUDGE_STORE_TRANSFORMS) l_store_transforms=trim(NUDGE_STORE_TRANSFORMS)=="YES" IMbkg=IM JMbkg=JM LMbkg=LM ! Validate grid ! ------------- call ESMF_GridValidate(GRID,RC=STATUS) VERIFY_(STATUS) call MAPL_GridGet(GRID, globalCellCountPerDim=DIMS, RC=STATUS) VERIFY_(STATUS) IMbkg_World=DIMS(1) JMbkg_World=DIMS(2) L_CUBE = JMbkg_World==6*IMbkg_World call ESMF_ClockGet(clock, currTime=currTime, rc=status) VERIFY_(STATUS) call ESMF_TimeGet(currTIME, timeString=DATE, RC=STATUS) VERIFY_(STATUS) call strToInt(DATE, nymd, nhms) call ESMF_CFIOstrTemplate ( AINCFILE, FILETMPL, 'GRADS', nymd=nymd, nhms=nhms, stat=STATUS ) VERIFY_(STATUS) if(MAPL_AM_I_ROOT()) then print * print *, 'Overwriting IAU-inc with ANA-Inc: ', trim(AINCFILE) endif call MAPL_GetResource( STATE, NX, Label="NX:", RC=status ) VERIFY_(STATUS) call MAPL_GetResource( STATE, NY, Label="NY:", RC=status ) VERIFY_(STATUS) call MAPL_GetResource( STATE, idum, 'ANALYZE_TS:', default=0, RC=STATUS ) VERIFY_(STATUS) IuseTS=idum/=0 if (.not.myANA%initialized) then call CFIO_Open ( AINCFILE, 1, fid, STATUS ) VERIFY_(STATUS) call CFIO_DimInquire ( fid, IMana_World, JMana_world, LM, nt, nvars, natts, STATUS ) VERIFY_(STATUS) call CFIO_Close ( fid, STATUS ) VERIFY_(STATUS) write(imstr,*) IMana_World write(jmstr,*) JMana_World if(JMana_World==6*IMana_World) then myANA%gridAnaName='PE'//trim(adjustl(imstr))//'x'//trim(adjustl(jmstr))//'-CF' myANA%GRIDana = AppGridCreateF(IMana_World,JMana_World,LM,NX,NY,STATUS) call WRITE_PARALLEL("Created cube myAna%GRIDana...") else myANA%gridAnaName='PC'//trim(adjustl(imstr))//'x'//trim(adjustl(jmstr))//'-DC' myANA%GRIDana = MAPL_LatLonGridCreate (Name=trim(myAna%gridAnaName), & NX = NX, & NY = NY, & IM_World = IMana_World, & JM_World = JMana_World, & LM_World = LM, & RC=STATUS ) call WRITE_PARALLEL("Created lat/lon myAna%GRIDana...") endif ! Validate grid ! ------------- call ESMF_GridValidate(myAna%GRIDana,RC=STATUS) VERIFY_(STATUS) call MAPL_GridGet(myANA%GRIDana, localCellCountPerDim=DIMS, RC=STATUS) VERIFY_(STATUS) myANA%IM = DIMS(1) myANA%JM = DIMS(2) myANA%LM = DIMS(3) ! Create transforms handles ! ------------------------- myANA%do_transforms = ( IMbkg_World /= IMana_World ) .or. & ( JMbkg_World /= JMana_World ) if (myANA%do_transforms) then if(MAPL_AM_I_ROOT()) then print * print *, 'Ana res: ', IMana_World,JMana_World,' Bkg res: ',IMbkg_World,JMbkg_World print * endif call MAPL_HorzTransformCreate( myAna%ANA2BKG, myAna%GRIDana, GRID , RC=STATUS ) VERIFY_(STATUS) call MAPL_HorzTransformCreate( myAna%BKG2ANA, GRID , myAna%GRIDana, RC=STATUS ) VERIFY_(STATUS) call WRITE_PARALLEL("Created transforms Ana2Bkg/Bkg2Ana ...") endif ! When remap required, need to get PHIS from model state onto analysis grid ! ------------------------------------------------------------------------- if (l_remap) then ! Get PHIS from background ! ------------------------ call ESMF_GridCompGet( GC, CONFIG = CF, RC=STATUS ) VERIFY_(STATUS) call ESMF_StateGet( GIM(SDYN), 'PHIS', FIELD, rc=STATUS ) VERIFY_(STATUS) Call GEOS_TopoGet ( CF, MEAN=FIELD, rc=STATUS ) VERIFY_(STATUS) call ESMF_StateGet( GIM(PHYS), 'PHIS', FIELD, rc=STATUS ) VERIFY_(STATUS) Call GEOS_TopoGet ( CF, MEAN=FIELD, rc=STATUS ) VERIFY_(STATUS) call ESMF_FieldGet (FIELD, localDE=0, farrayPtr=phis_bkg, rc = status) ! Derive PHIS as consequence of going from BKG to ANA back to BKG grid ! -------------------------------------------------------------------- allocate(myANA%phis_bkg(myANA%IM,myANA%JM),stat=STATUS);VERIFY_(STATUS) if (myANA%do_transforms) then allocate(qdum1(myANA%IM,myANA%JM,1)) allocate(qdum2(IMbkg,JMbkg,1)) qdum2(:,:,1)=phis_bkg call MAPL_HorzTransformRun (myANA%BKG2ANA, qdum2, qdum1, RC=STATUS );VERIFY_(STATUS) myANA%phis_bkg=qdum1(:,:,1) deallocate(qdum2) deallocate(qdum1) else myANA%phis_bkg=phis_bkg endif endif myAna%Initialized=.true. endif IMana=myANA%IM JMana=myANA%JM LMana=myANA%LM do_transforms = myANA%do_transforms if (do_transforms) then ANA2BKG => myANA%ANA2BKG BKG2ANA => myANA%BKG2ANA endif AINC_TIME(1) = nymd/10000 AINC_TIME(2) = mod(nymd,10000)/100 AINC_TIME(3) = mod(nymd,100) AINC_TIME(4) = nhms/10000 AINC_TIME(5) = mod(nhms,10000)/100 AINC_TIME(6) = mod(nhms,100) ! Set analysis increment time ! --------------------------- call ESMF_TimeSet( AincTime, YY = AINC_TIME(1), & MM = AINC_TIME(2), & DD = AINC_TIME(3), & H = AINC_TIME(4), & M = AINC_TIME(5), & S = AINC_TIME(6), rc=status ); VERIFY_(STATUS) ! Get MPI communicator ! -------------------- call ESMF_VMGetCurrent(vm, rc=status) VERIFY_(STATUS) call ESMF_VmGet(VM, mpicommunicator=vm_comm, rc=status) VERIFY_(STATUS) ! ***************************************************************************** ! **** READ Internal STATE (ie. ANA.ETA) from REPLAY File into BUNDLE **** ! ***************************************************************************** RBundle = ESMF_FieldBundleCreate( RC=STATUS) VERIFY_(STATUS) call ESMF_FieldBundleSet(RBundle, grid=myAna%GRIDana, rc=status) VERIFY_(STATUS) call MAPL_CFIORead ( AINCFILE, AincTime, Rbundle , RC=status) VERIFY_(STATUS) call ESMF_FieldBundleGet ( RBUNDLE, fieldCount=NQ, RC=STATUS ) VERIFY_(STATUS) ! Get pointer to background fields (current GCM state) ! ---------------------------------------------------- call MAPL_GetPointer( GEX(SDYN), AK,'AK', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), BK,'BK', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), delp_bkg,'DELP',RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), u_bkg,'U', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), v_bkg,'V', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), tv_bkg,'TV', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer( GEX(PHYS),o3_bkg,'O3PPMV', RC=STATUS) VERIFY_(STATUS) do K=1,NumFriendly NULLIFY(q_bkg) call ESMFL_BundleGetPointerToData(BUNDLE, Names(K), q_bkg, RC=STATUS ) VERIFY_(STATUS) if(Names(K)=='Q') exit enddo call MAPL_GetPointer( GEX(PHYS), ts_bkg, 'TS', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer( GEX(SDYN), ps_bkg, 'PS', RC=STATUS) VERIFY_(STATUS) ! Loop over GSI increment fields ! ------------------------------ if (l_nudge) then ! Bundle has analysis, leave it on its own grid allocate(gptr3d(IMana,JMana,LMana), stat=status );VERIFY_(STATUS) allocate(gptr2d(IMana,JMana), stat=status );VERIFY_(STATUS) allocate(dp_aux (IMana,JMana,LMana),stat=status );VERIFY_(STATUS) allocate(dts_aux (IMana,JMana), stat=status );VERIFY_(STATUS) allocate(du_aux (IMana,JMana,LMana),stat=status );VERIFY_(STATUS) allocate(dv_aux (IMana,JMana,LMana),stat=status );VERIFY_(STATUS) allocate(dt_aux (IMana,JMana,LMana),stat=status );VERIFY_(STATUS) allocate(dq_aux (IMana,JMana,LMana),stat=status );VERIFY_(STATUS) allocate(do3_aux (IMana,JMana,LMana),stat=status );VERIFY_(STATUS) allocate(dple_aux(IMana,JMana,0:LMana),stat=status );VERIFY_(STATUS) allocate(ple_bkg (IMbkg,JMbkg,0:LMbkg),stat=status );VERIFY_(STATUS) else ! Bundle has increment, therefore bring it to GCM grid allocate(gptr3d(IMbkg,JMbkg,LMbkg),stat=STATUS); VERIFY_(STATUS) allocate(gptr2d(IMbkg,JMbkg), stat=STATUS); VERIFY_(STATUS) allocate(qdum1 (IMbkg,JMbkg,1),stat=STATUS); VERIFY_(STATUS) allocate(qdum2 (IMana,JMana,1),stat=STATUS); VERIFY_(STATUS) endif allocate(phis_ana(IMana,JMana),stat=STATUS) VERIFY_(STATUS) allocate(dp_inc(size(gptr3d,1),size(gptr3d,2),size(gptr3d,3)),stat=STATUS) VERIFY_(STATUS) allocate(dps(size(gptr2d,1),size(gptr2d,2)),stat=STATUS) VERIFY_(STATUS) fromANA2BKG=(.not.l_nudge) .and. do_transforms do ni = 1, nq call ESMF_FieldBundleGet(RBUNDLE, ni, Field, __RC__ ) call ESMF_FieldGet(Field, NAME=NAME, dimCount = rank, __RC__ ) if (l_nudge) then if (.not.check_list_(NAME,ananames)) cycle else if (.not.check_list_(NAME,incnames)) cycle endif if (rank==2) then call ESMF_FieldGet(Field, farrayPtr=aptr2d, __RC__ ) if (fromANA2BKG) then qdum2(:,:,1)=aptr2d call MAPL_HorzTransformRun (ANA2BKG, qdum2, qdum1, RC=STATUS ) VERIFY_(STATUS) gptr2d=qdum1(:,:,1) else gptr2d=aptr2d endif if (fromANA2BKG) then if(trim(NAME)=='ts') dts_inc = gptr2d else if(trim(NAME)=='ts') dts_aux = gptr2d endif if(trim(NAME)=='ps') then dps=gptr2d endif if(trim(NAME)=='phis') then phis_ana=gptr2d endif else call ESMF_FieldGet(Field, farrayPtr=aptr3d, __RC__ ) if(trim(NAME)=='u'.or.trim(NAME)=='v') then if(trim(NAME)=='u') then allocate(uptr(IMana,JMana,LMana)) uptr=aptr3d cycle endif if(trim(NAME)=='v') then allocate(vptr(IMana,JMana,LMana)) vptr=aptr3d cycle endif endif if (fromANA2BKG) then call MAPL_HorzTransformRun (ANA2BKG, aptr3d, gptr3d, RC=STATUS ) VERIFY_(STATUS) else gptr3d=aptr3d endif if (fromANA2BKG) then if(trim(NAME)=='ozone') do3_inc=gptr3d if(trim(NAME)=='sphu' ) dq_inc =gptr3d if(trim(NAME)=='tv' ) dt_inc =gptr3d if(trim(NAME)=='delp' ) dp_inc =gptr3d else if(trim(NAME)=='ozone') do3_aux=gptr3d if(trim(NAME)=='sphu' ) dq_aux =gptr3d if(trim(NAME)=='tv' ) dt_aux =gptr3d if(trim(NAME)=='delp' ) dp_aux =gptr3d endif endif enddo if (fromANA2BKG) then deallocate(qdum2) deallocate(qdum1) endif ! U and V ! ------- ! could apply wind fix here ... but conversion to cubed ! TBD if (fromANA2BKG) then if( L_CUBE ) then call MAPL_HorzTransformRun (ANA2BKG, uptr, vptr, du_aux, dv_aux, RC=STATUS ) VERIFY_(STATUS) else call MAPL_HorzTransformRun (ANA2BKG, uptr, du_aux, RC=STATUS ) VERIFY_(STATUS) call MAPL_HorzTransformRun (ANA2BKG, vptr, dv_aux, RC=STATUS ) VERIFY_(STATUS) ! call POLEFIX ( du_aux,dv_aux,VM,GRID ) endif else du_aux=uptr dv_aux=vptr endif ! Calculate 3d-pressure change ! ----------------------------- if (.not.l_nudge) then ! increment dple_inc(:,:,0)=0.0 do L=1,LM dple_inc(:,:,L)=dp_inc(:,:,L) enddo endif ! Convert virtual temperature increment into dry temperature increment ! ------------------------------------------------------------------- EPS = MAPL_RVAP/MAPL_RGAS-1.0 if (.not.l_nudge) then ! in this case, using the background fields is not ! quite legitimate since these refer to the really ! current trajectory and not quite the original ! background used by the analysis dt_inc = dt_inc /(1.0+eps*q_bkg) - eps*dq_inc*tv_bkg/((1.0+eps*q_bkg)*(1.0+eps*q_bkg)) ! now dt is inc in dry temperature endif ! When nugding, RBundle carries full analysis state. ! In this case the follow will take place: ! 1. Convert present state of GCM to analysis grid ! 2. Calculate difference of analysis and present GCM state on analysis grid ! 3. Convert difference field back to GCM grid and overwrite analysis "tendencies" ! ----------------------------------------------------------------------------------- if (l_nudge) then allocate(u_ana (IMana,JMana, LMana),stat=STATUS); VERIFY_(STATUS) allocate(v_ana (IMana,JMana, LMana),stat=STATUS); VERIFY_(STATUS) allocate(tv_ana (IMana,JMana, LMana),stat=STATUS); VERIFY_(STATUS) allocate(q_ana (IMana,JMana, LMana),stat=STATUS); VERIFY_(STATUS) allocate(o3_ana (IMana,JMana, LMana),stat=STATUS); VERIFY_(STATUS) allocate(ple_ana(IMana,JMana,0:LMana),stat=STATUS); VERIFY_(STATUS) u_ana = du_aux v_ana = dv_aux q_ana = dq_aux tv_ana= dt_aux o3_ana= do3_aux ! Analyzed pressure edges ! ----------------------- ple_ana(:,:,0) = ak(0) do L=1,LM ple_ana(:,:,L) = ple_ana(:,:,L-1) + dp_aux(:,:,L) enddo ! Background pressure edges ! ------------------------- ple_bkg(:,:,0) = ak(0) do L=1,LM ple_bkg(:,:,L) = ple_bkg(:,:,L-1) + delp_bkg(:,:,L) enddo if (do_transforms) then ! Winds: if( L_CUBE ) then call MAPL_HorzTransformRun (BKG2ANA, u_bkg, v_bkg, du_aux, dv_aux, RC=STATUS ) VERIFY_(STATUS) else call MAPL_HorzTransformRun (BKG2ANA, u_bkg, du_aux, RC=STATUS ) VERIFY_(STATUS) call MAPL_HorzTransformRun (BKG2ANA, v_bkg, dv_aux, RC=STATUS ) VERIFY_(STATUS) ! call POLEFIX ( du_aux,dv_aux,VM,myANA%GRIDana ) endif ! Specific humidity: call MAPL_HorzTransformRun (BKG2ANA, q_bkg, dq_aux, RC=STATUS ) VERIFY_(STATUS) ! Pressure edges: call MAPL_HorzTransformRun (BKG2ANA, ple_bkg, dple_aux, RC=STATUS );VERIFY_(STATUS) ! Virtutal Temperature: call MAPL_HorzTransformRun (BKG2ANA, tv_bkg, dt_aux, RC=STATUS ) VERIFY_(STATUS) ! Ozone: call MAPL_HorzTransformRun (BKG2ANA, o3_bkg, do3_aux, RC=STATUS ) VERIFY_(STATUS) if (l_remap) then if(MAPL_AM_I_ROOT()) then print *, 'Remapping ANA to Internal State Topography' print * endif allocate(pk_ana (IMana,JMana, LMana),stat=STATUS);VERIFY_(STATUS) allocate(pke_ana(IMana,JMana,0:LMana),stat=STATUS);VERIFY_(STATUS) pke_ana(:,:,:) = ple_ana(:,:,:)**MAPL_KAPPA do L=1,lm pk_ana(:,:,L) = ( pke_ana(:,:,L)-pke_ana(:,:,L-1) ) & / ( MAPL_KAPPA*log(ple_ana(:,:,L)/ple_ana(:,:,L-1)) ) enddo allocate(thv_ana(IMana,JMana,LMana),stat=STATUS);VERIFY_(STATUS) thv_ana = tv_ana/pk_ana call myremap ( ple_ana, & u_ana, & v_ana, & thv_ana, & q_ana, & o3_ana, & phis_ana,myANA%phis_bkg,ak,bk,IMana,JMana,LMana ) ! Re-create ANA Dry Temperature ! ----------------------------- pke_ana(:,:,:) = ple_ana(:,:,:)**MAPL_KAPPA do L=1,LMana pk_ana(:,:,L) = ( pke_ana(:,:,L)-pke_ana(:,:,L-1) ) & / ( MAPL_KAPPA*log(ple_ana(:,:,L)/ple_ana(:,:,L-1)) ) enddo tv_ana= thv_ana*pk_ana deallocate(thv_ana,stat=STATUS);VERIFY_(STATUS) deallocate(pke_ana,stat=STATUS);VERIFY_(STATUS) deallocate(pk_ana ,stat=STATUS);VERIFY_(STATUS) endif ! ! Modify Vertically Integrated Mass-Divergence Increment ! ------------------------------------------------------ if ( l_windfix ) then if(MAPL_AM_I_ROOT()) then print *, 'Applying Mass Divergence Fix ...' print * endif method = 1 allocate ( vintdiva(IMana,JMana) ) allocate ( vintdivb(IMana,JMana) ) allocate ( vintdivc(IMana,JMana) ) call windfix ( u_ana,v_ana,ple_ana, & du_aux,dv_aux,dple_aux,IMana,JMana,LMana,VM,myANA%GRIDana,method, & vintdiva,vintdivb,vintdivc ) deallocate ( vintdivc ) deallocate ( vintdivb ) deallocate ( vintdiva ) endif ! Calculate wind increment on analysis grid du_aux = u_ana - du_aux dv_aux = v_ana - dv_aux ! Bring wind increments from analysis grid to model grid if( L_CUBE ) then call MAPL_HorzTransformRun (ANA2BKG, du_aux, dv_aux, du_inc, dv_inc, RC=STATUS ) VERIFY_(STATUS) else call MAPL_HorzTransformRun (ANA2BKG, du_aux, du_inc, RC=STATUS ) VERIFY_(STATUS) call MAPL_HorzTransformRun (ANA2BKG, dv_aux, dv_inc, RC=STATUS ) VERIFY_(STATUS) ! call POLEFIX ( du_inc,dv_inc,VM,myANA%GRIDana ) endif ! Calculate specific humdity increment on analysis grid dq_aux = q_ana - dq_aux ! Bring specific humidity increment from analysis grid to model grid call MAPL_HorzTransformRun (ANA2BKG, dq_aux, dq_inc, RC=STATUS ) ! Calculate increment to pressure thickness dple_aux = ple_ana - dple_aux call MAPL_HorzTransformRun (ANA2BKG, dple_aux, dple_inc, RC=STATUS );VERIFY_(STATUS) ! Calculate virtual temperature increment dt_aux = (tv_ana - dt_aux) ! virtual temperature increment call MAPL_HorzTransformRun (ANA2BKG, dt_aux, dt_inc, RC=STATUS ) ! Convert virtual temperature increment into dry temperature increment dt_inc = dt_inc/(1.0+eps*q_bkg) - eps*dq_inc*tv_bkg/((1.0+eps*q_bkg)*(1.0+eps*q_bkg)) ! dt_inc now has inc on dry temperature ! Calculate ozone increment on analysis grid do3_aux = o3_ana - do3_aux ! Bring specific humidity increment from analysis grid to model grid call MAPL_HorzTransformRun (ANA2BKG, do3_aux, do3_inc, RC=STATUS ) ! Calcuate T-skin increment if (IuseTS) then allocate(qdum2(IMana,JMana,1)) allocate(qdum1(IMbkg,JMbkg,1)) qdum1(:,:,1)=ts_bkg call MAPL_HorzTransformRun (BKG2ANA, qdum1, qdum2, RC=STATUS ) qdum2(:,:,1) = dts_aux - qdum2(:,:,1) call MAPL_HorzTransformRun (ANA2BKG, qdum2, qdum1, RC=STATUS ) VERIFY_(STATUS) dts_inc = qdum1(:,:,1) deallocate(qdum1) deallocate(qdum2) else dts_inc = 0.0 endif else ! both analysis and GCM are on the same grid (likely lat/lon) if(MAPL_AM_I_ROOT()) then print * print *, 'Handling same resolution case' print * endif ! Wind increment du_inc = du_aux - u_bkg dv_inc = dv_aux - v_bkg ! Specific humidity increment dq_inc = dq_aux - q_bkg ! Pressure increments dple_inc = ple_ana - ple_bkg ! Virtual Temperature dt_inc = (tv_ana - tv_bkg) ! virtual temperature increment dt_inc = dt_inc/(1.0+eps*q_bkg) - eps*dq_inc*tv_bkg/((1.0+eps*q_bkg)*(1.0+eps*q_bkg)) ! dry temperature increment ! Ozone do3_inc = do3_aux - o3_bkg ! Skin temperature if (IuseTS) then dts_inc = dts_aux - ts_bkg else dts_inc = 0.0 endif endif deallocate(ple_ana,stat=STATUS); VERIFY_(STATUS) deallocate(o3_ana, stat=STATUS); VERIFY_(STATUS) deallocate(q_ana, stat=STATUS); VERIFY_(STATUS) deallocate(tv_ana, stat=STATUS); VERIFY_(STATUS) deallocate(v_ana, stat=STATUS); VERIFY_(STATUS) deallocate(u_ana, stat=STATUS); VERIFY_(STATUS) endif ! ! Clean up ! -------- call MAPL_FieldBundleDestroy ( RBundle, RC=STATUS) VERIFY_(STATUS) if(l_nudge) then deallocate(ple_bkg, stat=STATUS); VERIFY_(STATUS) deallocate(dts_aux ,stat=STATUS); VERIFY_(STATUS) deallocate(du_aux ,stat=STATUS); VERIFY_(STATUS) deallocate(dv_aux ,stat=STATUS); VERIFY_(STATUS) deallocate(dp_aux ,stat=STATUS); VERIFY_(STATUS) deallocate(dt_aux ,stat=STATUS); VERIFY_(STATUS) deallocate(dq_aux ,stat=STATUS); VERIFY_(STATUS) deallocate(do3_aux ,stat=STATUS); VERIFY_(STATUS) deallocate(dple_aux,stat=STATUS); VERIFY_(STATUS) endif if(associated(phis_ana))deallocate(phis_ana) if(associated(gptr2d))deallocate(gptr2d) if(associated(gptr3d))deallocate(gptr3d) if(allocated(dps)) deallocate(dps) if(allocated(dp_inc))deallocate(dp_inc) if(allocated(uptr)) deallocate(uptr) if(allocated(vptr)) deallocate(vptr) if ( (.not.l_store_transforms) ) then if (associated(myANA%phis_bkg)) then deallocate(myANA%phis_bkg,stat=STATUS);VERIFY_(STATUS) nullify(myANA%phis_bkg) endif if( do_transforms ) then call MAPL_HorzTransformDestroy(BKG2ANA, RC=STATUS);VERIFY_(STATUS) call MAPL_HorzTransformDestroy(ANA2BKG, RC=STATUS);VERIFY_(STATUS) call WRITE_PARALLEL("Destroyed transforms ANA2BKG/BKG2ANA") endif call ESMF_GridDestroy(myAna%GRIDana, rc=status);VERIFY_(STATUS) call WRITE_PARALLEL("Destroyed myAna%GRIDana") myANA%do_transforms=.false. myANA%initialized=.false. endif if( LAST_CORRECTOR .and. myANA%initialized ) then if (associated(myANA%phis_bkg)) then deallocate(myANA%phis_bkg,stat=STATUS);VERIFY_(STATUS) nullify(myANA%phis_bkg) endif if( do_transforms ) then call MAPL_HorzTransformDestroy(BKG2ANA, RC=STATUS);VERIFY_(STATUS) call MAPL_HorzTransformDestroy(ANA2BKG, RC=STATUS);VERIFY_(STATUS) call WRITE_PARALLEL("Destroyed transforms ANA2BKG/BKG2ANA") endif call ESMF_GridDestroy(myAna%GRIDana, rc=status);VERIFY_(STATUS) call WRITE_PARALLEL("Destroyed myAna%GRIDana") myANA%do_transforms=.false. myANA%initialized=.false. endif RETURN_(ESMF_SUCCESS) end subroutine update_ainc_ logical function check_list_(name,vars) implicit none character(len=*) :: name character(len=*) :: vars(:) integer ii check_list_=.false. do ii = 1,size(vars) if(trim(name)==trim(vars(ii))) then check_list_=.true. exit endif enddo end function check_list_ subroutine dfi_coeffs (FilterOpt,DT,TC,TAUANL,nsteps,dfi) ! This subroutine belongs to GEOS_Shared, but for now it lives here implicit none character(len=*),intent(in):: FilterOpt ! IAU or DFI real, intent(in) :: DT ! model time step real, intent(in) :: TC ! cutoff time (typically 21600 sec) real, intent(in) :: TAUANL ! IAU coefficient integer,intent(in) :: nsteps ! number of steps TC/DT+1 real, intent(out) :: dfi(nsteps) real pi,arg,wc integer n,i,k,nhlf,np1 ! If regular IAU, simply set to constant and return ! ------------------------------------------------- if (trim(FilterOpt)=='IAU' .or. trim(FilterOpt)=='iau' ) then dfi=1.0/TAUANL return endif ! Calculate DFI coefficients ! -------------------------- if (abs(TAUANL-TC) < 1.e-5) then if(MAPL_AM_I_ROOT()) then print *, 'Warning: Inconsistent IAU/filtering function' print *, 'Warning: Check TAUANL and POFFSET' print *, 'Warning: Ignoring TAUANL ...' endif endif pi = 4.0*atan(1.) nhlf = (nsteps+1)/2 do k = 1, nhlf-1 n = k-nhlf arg = n*pi/nhlf wc = sin(arg)/arg ! Lanczos window dfi(k) = wc*sin(n*2.0*pi*DT/TC)/(n*pi) end do dfi(nhlf) = 2*DT/TC do i = nhlf+1, nsteps dfi(i) = dfi(nsteps-i+1) end do ! Normalize coefficients ! ---------------------- dfi = dfi/sum(dfi) dfi = dfi/DT ! remember: dynamics multiplies by DT end subroutine dfi_coeffs end subroutine Run function my_nearest_time(TimeLeft, TimeRight, TimeNow, RC) result(TimeNearest) type(ESMF_Time), intent(IN) :: TimeLeft type(ESMF_Time), intent(IN) :: TimeRight type(ESMF_Time), intent(IN) :: TimeNow integer, optional, intent(OUT) :: RC type(ESMF_Time) :: TimeNearest if (TimeLeft==TimeRight) then TimeNearest=TimeLeft if(present(RC)) then RC = 0 endif return endif if (TimeNow<=TimeLeft) then TimeNearest=TimeLeft if(present(RC)) then RC = 0 endif return endif if (TimeNow>=TimeRight) then TimeNearest=TimeRight if(present(RC)) then RC = 0 endif return endif if ((TimeNow-TimeLeft) < (TimeNow-TimeRight)) then TimeNearest=TimeLeft else TimeNearest=TimeRight endif if(present(RC)) then RC = 0 endif end function my_nearest_time end module GEOS_AgcmGridCompMod