! +-======-+ ! 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: GEOS_mkiauGridComp.F90,v 1.17 2012-11-19 21:08:23 ltakacs Exp $ #include "MAPL_Generic.h" !============================================================================= !BOP ! !MODULE: GEOS_mkiau -- A Module to compute the IAU forcing ! !INTERFACE: module GEOS_mkiauGridCompMod ! !USES: use ESMF_Mod use MAPL_Mod use G3_MPI_Util_Mod use GEOS_RemapMod, only: myremap => remap use m_mpif90, only: mp_comm_world use m_chars, only: uppercase implicit none private ! !PUBLIC MEMBER FUNCTIONS: public SetServices !============================================================================= ! !DESCRIPTION: ! ! !EOP 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: This version uses the MAPL_GenericSetServices. This function sets ! the Initialize and Finalize services, as well as allocating ! our instance of a generic state and putting it in the ! gridded component (GC). Here we only need to set the run method and ! add the state variable specifications (also generic) to our instance ! of the generic state. This is the way our true state variables get into ! the ESMF_State INTERNAL, which is in the MAPL_MetaComp. !EOP !============================================================================= ! ! ErrLog Variables character(len=ESMF_MAXSTR) :: IAm integer :: STATUS integer :: OFFLINE character(len=ESMF_MAXSTR) :: COMP_NAME type (MAPL_MetaComp), pointer :: MAPL !============================================================================= ! Begin... ! Get my name and set-up traceback handle ! --------------------------------------- Iam = 'SetServices' call ESMF_GridCompGet( GC, NAME=COMP_NAME, RC=STATUS ) VERIFY_(STATUS) Iam = trim(COMP_NAME) // Iam ! Retrieve the pointer to the state !---------------------------------- call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS) VERIFY_(STATUS) ! Set the Run entry point ! ----------------------- call MAPL_GridCompSetEntryPoint ( gc, ESMF_SETRUN, Run, & RC=STATUS) VERIFY_(STATUS) ! Set the state variable specs. ! ----------------------------- call MAPL_GetResource(MAPL, OFFLINE, LABEL="REPLAY_OFFLINE:", default=0, RC=status) VERIFY_(STATUS) ! !IMPORT STATE: ! -------------- call MAPL_AddImportSpec ( GC, & SHORT_NAME = 'PHIS', & LONG_NAME = 'surface geopotential height', & UNITS = 'm2 sec-2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( GC, & SHORT_NAME = 'AK', & LONG_NAME = 'hybrid_sigma_pressure_a', & UNITS = 'Pa', & DIMS = MAPL_DimsVertOnly, & VLOCATION = MAPL_VLocationEdge, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( GC, & SHORT_NAME = 'BK', & LONG_NAME = 'hybrid_sigma_pressure_b', & UNITS = 'Pa', & DIMS = MAPL_DimsVertOnly, & VLOCATION = MAPL_VLocationEdge, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'PLE', & LONG_NAME = 'air_pressure', & UNITS = 'Pa', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationEdge, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'T', & LONG_NAME = 'air_temperature', & UNITS = 'K', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'U', & LONG_NAME = 'eastward_wind', & UNITS = 'm s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'V', & LONG_NAME = 'northward_wind', & UNITS = 'm s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'O3PPMV', & LONG_NAME = 'ozone_volume_mixing_ratio', & UNITS = 'ppmv', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec ( gc, & SHORT_NAME = 'TS', & LONG_NAME = 'skin_temperature', & UNITS = 'K', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) ! we need the analysis tracers' bundle (for QV and O3) ! ---------------------------------------------------- if ( OFFLINE==0 ) then call MAPL_AddImportSpec(GC, & SHORT_NAME = 'TRANA', & LONG_NAME = 'analyzed_quantities', & UNITS = 'X', & DATATYPE = MAPL_BundleItem, & RC=STATUS ) VERIFY_(STATUS) else call MAPL_AddImportSpec(GC, & SHORT_NAME = 'QV', & LONG_NAME = 'water_vapor_specific_humdity', & UNITS = 'kg/kg', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) endif ! !EXPORT STATE: ! -------------- call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DUDT', & LONG_NAME = 'eastward_wind_analysis_increment', & UNITS = 'm s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DVDT', & LONG_NAME = 'northward_wind_analysis_increment', & UNITS = 'm s-1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DTDT', & LONG_NAME = 'temperature_analysis_increment', & UNITS = 'K', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DPEDT', & LONG_NAME = 'edge_pressure_analysis_increment', & UNITS = 'Pa', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationEdge, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( 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_AddExportSpec ( gc, & SHORT_NAME = 'DO3DT', & LONG_NAME = 'ozone_analysis_increment', & UNITS = 'ppmv', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec ( gc, & SHORT_NAME = 'DTSDT', & LONG_NAME = 'skin_temparature_increment', & UNITS = 'K', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, RC=STATUS ) VERIFY_(STATUS) ! Internal State (None) ! --------------------- ! Set the Profiling timers ! ------------------------ call MAPL_TimerAdd(GC, name="INITIALIZE" ,RC=STATUS) VERIFY_(STATUS) call MAPL_TimerAdd(GC, name="DRIVER" ,RC=STATUS) VERIFY_(STATUS) call MAPL_TimerAdd(GC, name="-INTR" ,RC=STATUS) VERIFY_(STATUS) ! Set generic init and final methods ! ---------------------------------- call MAPL_GenericSetServices ( gc, RC=STATUS) VERIFY_(STATUS) RETURN_(ESMF_SUCCESS) end subroutine SetServices !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !BOP ! ! IROUTINE: RUN -- Run method for the MAKEIAU 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: type ( dynamics_lattice_type ) lattice ! ! DESCRIPTION: This version uses the MAPL_GenericSetServices. This function sets ! the Initialize and Finalize services, as well as allocating !EOP ! ErrLog Variables character(len=ESMF_MAXSTR) :: IAm integer :: STATUS character(len=ESMF_MAXSTR) :: COMP_NAME ! Local derived type aliases type (MAPL_MetaComp), pointer :: MAPL type (ESMF_State ) :: INTERNAL integer :: IM, JM, LM real, pointer, dimension(:,:,:) :: u_bkg, du real, pointer, dimension(:,:,:) :: v_bkg, dv real, pointer, dimension(:,:,:) :: t_bkg, dt real, pointer, dimension(:,:,:) :: q_bkg, dq real, pointer, dimension(:,:,:) :: o3_bkg, do3 real, pointer, dimension(:,:,:) :: ple_bkg, dple real, pointer, dimension(:,:) :: ts_bkg, dts real, pointer, dimension(:,:) ::phis_bkg real, pointer, dimension(:) :: ak,bk real, pointer, dimension(:,:) ::phis_ana real, pointer, dimension(:,:) :: ts_ana real, pointer, dimension(:,:) :: ps_ana real, pointer, dimension(:,:,:) :: dp_ana real, pointer, dimension(:,:,:) :: u_ana real, pointer, dimension(:,:,:) :: v_ana real, pointer, dimension(:,:,:) :: t_ana real, pointer, dimension(:,:,:) :: thv_ana real, pointer, dimension(:,:,:) :: q_ana real, pointer, dimension(:,:,:) :: o3_ana real, allocatable, dimension(:,:,:) :: ple_ana real, allocatable, dimension(:,:,:) :: pk_ana real, allocatable, dimension(:,:,:) :: pke_ana real, allocatable, dimension(:,:) :: pl, alf real, parameter :: EPS = MAPL_RVAP/MAPL_RGAS-1.0 character(len=ESMF_MAXSTR) :: FILETMPL character(len=ESMF_MAXSTR) :: cremap type(ESMF_FieldBundle) :: RBUNDLE type(ESMF_Field) :: FIELD integer :: K,KTS,NQ character(len=ESMF_MAXSTR), ALLOCATABLE :: RNAMES(:) character(len=ESMF_MAXSTR) :: STRING real, pointer :: temp3d(:,:,:) real, pointer :: temp2d(:,:) type(ESMF_FieldBundle) :: bundle type(ESMF_Grid) :: grid type(ESMF_Time) :: currtime real :: DAMPBEG, DAMPEND integer :: L integer :: NX,NY,IMG,JMG integer :: method integer :: OFFLINE logical :: doremap logical :: FOUND logical :: ANALYZE_TS logical :: REPLAY_U, REPLAY_V, REPLAY_T, REPLAY_QV, REPLAY_TS, REPLAY_O3, REPLAY_P character(len=ESMF_MAXSTR) :: ALIAS_U, ALIAS_V, ALIAS_T, ALIAS_QV, ALIAS_TS, ALIAS_O3, ALIAS_PS, ALIAS_DP character(len=ESMF_MAXSTR) :: ALIAS_TH, ALIAS_TV, ALIAS_THV character(len=ESMF_MAXSTR) :: ALIAS_PHIS logical :: first data first /.true./ !============================================================================= ! 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) // Iam ! Retrieve the pointer to the state !---------------------------------- call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS) VERIFY_(STATUS) ! Local aliases to the state, grid, and configuration ! --------------------------------------------------- call MAPL_TimerOn(MAPL,"TOTAL") ! Get parameters from generic state. !----------------------------------- call MAPL_Get( MAPL, IM=IM, JM=JM, LM=LM, & INTERNAL_ESMF_STATE=INTERNAL, & RC=STATUS ) VERIFY_(STATUS) call ESMF_GridCompGet(GC, grid=grid, rc=status) VERIFY_(STATUS) ! Get Resource Parameters !------------------------ call MAPL_GetResource(MAPL, FILETMPL, LABEL="REPLAY_FILE:", RC=status) VERIFY_(STATUS) call MAPL_GetResource(MAPL, CREMAP, LABEL="REPLAY_REMAP:", default="yes", RC=status) VERIFY_(STATUS) call MAPL_GetResource(MAPL, OFFLINE, LABEL="REPLAY_OFFLINE:", default=0, RC=status) VERIFY_(STATUS) call MAPL_GetResource(MAPL, KTS, Label="ANALYZE_TS:", default=0, RC=STATUS) VERIFY_(STATUS) ANALYZE_TS = (KTS /= 0) call MAPL_GetResource(MAPL, K, Label="REPLAY_TS:", default=kts, RC=STATUS) VERIFY_(STATUS) REPLAY_TS = (K == 1) call MAPL_GetResource(MAPL, K, Label="REPLAY_P:", default=1, RC=STATUS) VERIFY_(STATUS) REPLAY_P = (K == 1) call MAPL_GetResource(MAPL, K, Label="REPLAY_U:", default=1, RC=STATUS) VERIFY_(STATUS) REPLAY_U = (K == 1) call MAPL_GetResource(MAPL, K, Label="REPLAY_V:", default=1, RC=STATUS) VERIFY_(STATUS) REPLAY_V = (K == 1) call MAPL_GetResource(MAPL, K, Label="REPLAY_T:", default=1, RC=STATUS) VERIFY_(STATUS) REPLAY_T = (K == 1) call MAPL_GetResource(MAPL, K, Label="REPLAY_QV:", default=1, RC=STATUS) VERIFY_(STATUS) REPLAY_QV = (K == 1) call MAPL_GetResource(MAPL, K, Label="REPLAY_O3:", default=1, RC=STATUS) VERIFY_(STATUS) REPLAY_O3 = (K == 1) call MAPL_GetResource(MAPL, ALIAS_PHIS, Label="REPLAY_PHIS_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, ALIAS_TS, Label="REPLAY_TS_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, ALIAS_PS, Label="REPLAY_PS_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, ALIAS_DP, Label="REPLAY_DP_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, ALIAS_U, Label="REPLAY_U_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, ALIAS_V, Label="REPLAY_V_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, ALIAS_T, Label="REPLAY_T_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, ALIAS_QV, Label="REPLAY_QV_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, ALIAS_O3, Label="REPLAY_O3_ALIAS:", default="NULL", RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, DAMPBEG, LABEL="REPLAY_DAMPBEG:", default=20.0, RC=status) VERIFY_(STATUS) call MAPL_GetResource(MAPL, DAMPEND, LABEL="REPLAY_DAMPEND:", default=40.0, RC=status) VERIFY_(STATUS) ASSERT_(DAMPBEG.le.DAMPEND ) CREMAP = uppercase(CREMAP) ! ********************************************************************** ! **** READ Internal STATE from REPLAY File **** ! ********************************************************************** call ESMF_ClockGet(clock, currTime=currTime, rc=status) VERIFY_(STATUS) RBundle = ESMF_FieldBundleCreate( GRID, RC=STATUS) VERIFY_(STATUS) call MAPL_CFIORead ( FILETMPL, currTime, Rbundle , RC=status) VERIFY_(STATUS) call ESMF_FieldBundleGet ( RBUNDLE, fieldCount=NQ, RC=STATUS ) VERIFY_(STATUS) if( .not.allocated( rnames ) ) then allocate( RNAMES(NQ),STAT=STATUS ) VERIFY_(STATUS) call ESMF_FieldBundleGet ( RBUNDLE, nameList=RNAMES, rc=STATUS ) VERIFY_(STATUS) if( first ) then write(STRING,'(A,I)') "REPLAY File Variables, NQ: ", nq call WRITE_PARALLEL( trim(STRING) ) do k=1,nq call WRITE_PARALLEL( trim(RNAMES(K)) ) enddo first = .false. endif endif ! ********************************************************************** ! **** Get Pointers to BKG Import Data **** ! ********************************************************************** call MAPL_GetPointer(import, u_bkg, 'U', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, v_bkg, 'V', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, t_bkg, 'T', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, ple_bkg, 'PLE', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, o3_bkg, 'O3PPMV', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import,phis_bkg, 'PHIS', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, ak, 'AK', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, bk, 'BK', RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(import, ts_bkg, 'TS', RC=STATUS) VERIFY_(STATUS) if ( OFFLINE==0 ) then call ESMF_StateGet(import, 'TRANA', bundle, rc=status) VERIFY_(STATUS) call ESMFL_BundleGetPointerToData(bundle, 'Q', q_bkg, RC=STATUS) VERIFY_(STATUS) else call MAPL_GetPointer(import, q_bkg, 'QV', RC=STATUS) VERIFY_(STATUS) endif ! ********************************************************************** ! **** Get Pointers to ANA Replay Data **** ! ********************************************************************** allocate ( phis_ana(IM,JM) ) allocate ( ts_ana(IM,JM) ) allocate ( ps_ana(IM,JM) ) allocate ( dp_ana(IM,JM, LM) ) allocate ( u_ana(IM,JM, LM) ) allocate ( v_ana(IM,JM, LM) ) allocate ( t_ana(IM,JM, LM) ) allocate ( thv_ana(IM,JM, LM) ) allocate ( q_ana(IM,JM, LM) ) allocate ( o3_ana(IM,JM, LM) ) allocate ( pk_ana(IM,JM, LM) ) allocate ( ple_ana(IM,JM,0:LM) ) allocate ( pke_ana(IM,JM,0:LM) ) doremap = trim(cremap).eq.'YES' ! Surface Temperature !-------------------- if( REPLAY_TS ) then FOUND = .false. do k=1,nq if( match('ts',alias_ts,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP2D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then ts_ana = TEMP2D FOUND = .true. exit endif endif enddo if( .not.FOUND ) then write(STRING,'(A)') "REPLAY Variable: TS Not Found!" call WRITE_PARALLEL( trim(STRING) ) RETURN_(ESMF_FAILURE) endif else ts_ana = ts_bkg endif ! Surface Geopotential Heights !----------------------------- FOUND = .false. do k=1,nq if( match('phis',alias_phis,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP2D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then phis_ana = TEMP2D FOUND = .true. exit endif endif enddo if( .not.FOUND ) then phis_ana = phis_bkg doremap =.false. endif ! Surface Pressure !----------------- if( REPLAY_P ) then FOUND = .false. do k=1,nq if( match('ps',alias_ps,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP2D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then ps_ana = TEMP2D FOUND = .true. exit endif endif enddo if( .not.FOUND ) then write(STRING,'(A)') "REPLAY Variable: PS Not Found!" call WRITE_PARALLEL( trim(STRING) ) RETURN_(ESMF_FAILURE) endif else ps_ana = ple_bkg(:,:,LM) endif ! Pressure Thickness !------------------- if( REPLAY_P ) then FOUND = .false. do k=1,nq if( match('dp',alias_dp,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then dp_ana = TEMP3D FOUND = .true. exit endif endif enddo if( .not.FOUND ) then write(STRING,'(A)') "REPLAY Variable: DP Not Found!" call WRITE_PARALLEL( trim(STRING) ) RETURN_(ESMF_FAILURE) endif else do L=1,lm dp_ana(:,:,L) = ple_bkg(:,:,L)-ple_bkg(:,:,L-1) end do endif ! U-Wind !------- if( REPLAY_U ) then FOUND = .false. do k=1,nq if( match('u',alias_u,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then u_ana = TEMP3D FOUND = .true. exit endif endif enddo if( .not.FOUND ) then write(STRING,'(A)') "REPLAY Variable: U Not Found!" call WRITE_PARALLEL( trim(STRING) ) RETURN_(ESMF_FAILURE) endif else u_ana = u_bkg endif ! V-Wind !------- if( REPLAY_V ) then FOUND = .false. do k=1,nq if( match('v',alias_v,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then v_ana = TEMP3D FOUND = .true. exit endif endif enddo if( .not.FOUND ) then write(STRING,'(A)') "REPLAY Variable: V Not Found!" call WRITE_PARALLEL( trim(STRING) ) RETURN_(ESMF_FAILURE) endif else v_ana = v_bkg endif ! Moisture Variable !------------------ if( REPLAY_QV ) then FOUND = .false. do k=1,nq if( match('qv',alias_qv,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then q_ana = TEMP3D FOUND = .true. exit endif endif enddo if( .not.FOUND ) then write(STRING,'(A)') "REPLAY Variable: QV Not Found!" call WRITE_PARALLEL( trim(STRING) ) RETURN_(ESMF_FAILURE) endif else q_ana = q_bkg endif ! Ozone !------ if( REPLAY_O3 ) then FOUND = .false. do k=1,nq if( match('o3',alias_o3,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then o3_ana = TEMP3D FOUND = .true. exit endif endif enddo if( .not.FOUND ) then write(STRING,'(A)') "REPLAY Variable: O3 Not Found!" call WRITE_PARALLEL( trim(STRING) ) RETURN_(ESMF_FAILURE) endif else o3_ana = o3_bkg endif ! Construct Pressure Variables ! ---------------------------- ple_ana(:,:,0) = ak(0) + ps_ana(:,:)*bk(0) do L=1,lm-1 ple_ana(:,:,L) = ple_ana(:,:,L-1) + dp_ana(:,:,L) end do ple_ana(:,:,LM) = ps_ana 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 ! Temperature Variable !--------------------- if( REPLAY_T ) then FOUND = .false. do k=1,nq if( match('t',alias_t,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then t_ana = TEMP3D FOUND = .true. exit endif endif if( match('tv',alias_tv,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then t_ana = TEMP3D/(1.0+eps*q_ana) FOUND = .true. exit endif endif if( match('th',alias_th,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then t_ana = TEMP3D*pk_ana FOUND = .true. exit endif endif if( match('thv',alias_thv,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then t_ana = TEMP3D*pk_ana/(1.0+eps*q_ana) FOUND = .true. exit endif endif enddo if( .not.FOUND ) then write(STRING,'(A)') "REPLAY Variable: T Not Found!" call WRITE_PARALLEL( trim(STRING) ) RETURN_(ESMF_FAILURE) endif else t_ana = t_bkg endif ! Test for Re-Mapping ! ------------------- if (doremap ) then if( count( phis_ana.ne.phis_bkg ).ne.0 ) then if(MAPL_AM_I_ROOT()) then print * print *, 'Remapping ANA Data to BKG Topography' print * endif ! Create ANA Virtual Potential Temperature ! ---------------------------------------- thv_ana = t_ana*(1.0+eps*q_ana)/pk_ana call myremap ( ple_ana, & u_ana, & v_ana, & thv_ana, & q_ana, & o3_ana, & phis_ana,phis_bkg,ak,bk,im,jm,lm ) ! Create ANA Dry Temperature ! -------------------------- 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 t_ana = thv_ana*pk_ana/(1.0+eps*q_ana) else if(MAPL_AM_I_ROOT()) then print * print *, 'ANA and BKG Topography match, No Remapping Necessary' print * endif endif else if(MAPL_AM_I_ROOT()) then print * print *, 'Not Remapping ANA Data to BKG Topography' print * endif endif ! ********************************************************************** ! **** Blend ANA and BKG Variables between DAMPBEG and DAMPEND **** ! ********************************************************************** if( DAMPBEG.ne.DAMPEND ) then call blend ( ple_ana,u_ana,v_ana,t_ana,q_ana,o3_ana, & ple_bkg,u_bkg,v_bkg,t_bkg,q_bkg,o3_bkg, im,jm,lm, DAMPBEG,DAMPEND ) call MAPL_GetResource(MAPL, NX, Label="NX:", RC=status ) VERIFY_(STATUS) call MAPL_GetResource(MAPL, NY, Label="NY:", RC=status ) VERIFY_(STATUS) call MAPL_GetResource(MAPL, IMG, Label="AGCM_IM:",RC=status ) VERIFY_(STATUS) call MAPL_GetResource(MAPL, JMG, Label="AGCM_JM:",RC=status ) VERIFY_(STATUS) ! Modify Vertically Integrated Mass-Divergence Increment (FV-Lat/Lon Core & IMG ! <= 576) ! ------------------------------------------------------------------------------------- if( JMG.ne.6*IMG .and. IMG.le.576 ) then call create_dynamics_lattice ( lattice,mp_comm_world,IMG,JMG,lm,NX,NY) method = 1 call windfix ( u_ana,v_ana,ple_ana, & u_bkg,v_bkg,ple_bkg,im,jm,lm,lattice,'A',method ) call destroy_dynamics_lattice (lattice) endif endif ! ********************************************************************** ! **** Create A-Grid IAU Increment **** ! ********************************************************************** call MAPL_GetPointer(export, du, 'DUDT', alloc=.TRUE., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(export, dv, 'DVDT', alloc=.TRUE., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(export, dt, 'DTDT', alloc=.TRUE., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(export, dq, 'DQVDT', alloc=.TRUE., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(export, do3, 'DO3DT', alloc=.TRUE., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(export, dple, 'DPEDT', alloc=.TRUE., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(export, dts, 'DTSDT', alloc=.TRUE., RC=STATUS) VERIFY_(STATUS) if( REPLAY_U ) then du = u_ana-u_bkg else du = 0.0 endif if( REPLAY_V ) then dv = v_ana-v_bkg else dv = 0.0 endif if( REPLAY_T ) then dt = t_ana-t_bkg else dt = 0.0 endif if( REPLAY_QV ) then dq = q_ana-q_bkg else dq = 0.0 endif if( REPLAY_O3 ) then do3 = o3_ana-o3_bkg else do3 = 0.0 endif if( REPLAY_P ) then dple = ple_ana-ple_bkg else dple = 0.0 endif if( REPLAY_TS ) then dts = ts_ana-ts_bkg else dts = 0.0 endif ! Clean-Up ! -------- do k=1,nq call ESMF_FieldBundleGet ( RBundle, k, Field, RC=STATUS) VERIFY_(STATUS) call MAPL_FieldF90Deallocate ( Field, RC=STATUS) VERIFY_(STATUS) call ESMF_FieldDestroy ( Field, RC=STATUS) VERIFY_(STATUS) enddo call ESMF_FieldBundleDestroy ( RBundle, RC=STATUS) VERIFY_(STATUS) deallocate ( RNAMES ) deallocate ( ple_ana ) deallocate ( phis_ana ) deallocate ( ts_ana ) deallocate ( u_ana ) deallocate ( v_ana ) deallocate ( t_ana ) deallocate ( q_ana ) deallocate ( o3_ana ) deallocate ( ps_ana ) deallocate ( dp_ana ) deallocate ( pk_ana ) deallocate ( pke_ana ) deallocate ( thv_ana ) call MAPL_TimerOff(MAPL,"TOTAL") RETURN_(ESMF_SUCCESS) end subroutine RUN subroutine blend ( plea,ua,va,ta,qa,oa, & pleb,ub,vb,tb,qb,ob, im,jm,lm, pabove,pbelow ) implicit none integer im,jm,lm real plea(im,jm,lm+1) real ua(im,jm,lm) real va(im,jm,lm) real ta(im,jm,lm) real qa(im,jm,lm) real oa(im,jm,lm) real pleb(im,jm,lm+1) real ub(im,jm,lm) real vb(im,jm,lm) real tb(im,jm,lm) real qb(im,jm,lm) real ob(im,jm,lm) real pabove,pbelow ! Locals ! ------ real pkea(im,jm,lm+1) real pkeb(im,jm,lm+1) real phia(im,jm,lm+1) real phib(im,jm,lm+1) real thva(im,jm,lm) real thvb(im,jm,lm) real pka(im,jm,lm) real pkb(im,jm,lm) real alf,eps,p integer i,j,L eps = MAPL_RVAP/MAPL_RGAS-1.0 pkea = plea**MAPL_KAPPA pkeb = pleb**MAPL_KAPPA do L=1,lm pka(:,:,L) = ( pkea(:,:,L+1)-pkea(:,:,L) ) / ( MAPL_KAPPA*log(plea(:,:,L+1)/plea(:,:,L)) ) pkb(:,:,L) = ( pkeb(:,:,L+1)-pkeb(:,:,L) ) / ( MAPL_KAPPA*log(pleb(:,:,L+1)/pleb(:,:,L)) ) thva(:,:,L) = ta(:,:,L)*(1.0+eps*qa(:,:,L)) / pka(:,:,L) thvb(:,:,L) = tb(:,:,L)*(1.0+eps*qb(:,:,L)) / pkb(:,:,L) enddo phia(:,:,lm+1) = 0.0 phib(:,:,lm+1) = 0.0 do L=lm,1,-1 phia(:,:,L) = phia(:,:,L+1) + MAPL_CP*thva(:,:,L)*( pkea(:,:,L+1)-pkea(:,:,L) ) phib(:,:,L) = phib(:,:,L+1) + MAPL_CP*thvb(:,:,L)*( pkeb(:,:,L+1)-pkeb(:,:,L) ) enddo ! Blend mid-level u,v,q and o3 ! ---------------------------- do L=1,lm do j=1,jm do i=1,im p = 0.5*( plea(i,j,L)+plea(i,j,L+1) ) if( p.le.pabove ) then alf = 0.0 else if( p.gt.pabove .and. p.le.pbelow ) then alf = (p-pabove)/(pbelow-pabove) else alf = 1.0 endif ua(i,j,L) = ub(i,j,L) + alf*( ua(i,j,L)- ub(i,j,L) ) va(i,j,L) = vb(i,j,L) + alf*( va(i,j,L)- vb(i,j,L) ) qa(i,j,L) = qb(i,j,L) + alf*( qa(i,j,L)- qb(i,j,L) ) oa(i,j,L) = ob(i,j,L) + alf*( oa(i,j,L)- ob(i,j,L) ) enddo enddo enddo ! Blend edge-level phi ! -------------------- do L=1,lm+1 do j=1,jm do i=1,im p = plea(i,j,L) if( p.le.pabove ) then alf = 0.0 else if( p.gt.pabove .and. p.le.pbelow ) then alf = (p-pabove)/(pbelow-pabove) else alf = 1.0 endif phia(i,j,L) = phib(i,j,L) + alf*( phia(i,j,L)-phib(i,j,L) ) enddo enddo enddo ! Compute T based on blended phi ! ------------------------------ do L=1,lm ta(:,:,L) = ( phia(:,:,L)-phia(:,:,L+1) )/( pkea(:,:,L+1)-pkea(:,:,L) ) & / (MAPL_CP*(1.0+eps*qa(:,:,L))) * pka(:,:,L) enddo return end subroutine blend function match (replay_name,replay_alias,replay_var) character(*) :: replay_name,replay_alias,replay_var character(len=ESMF_MAXSTR) :: name,alias,var logical match match = .false. name = uppercase( trim(replay_name ) ) alias = uppercase( trim(replay_alias) ) var = uppercase( trim(replay_var ) ) if( trim(var) == trim(alias) ) match = .true. if( trim(name) == 'U' ) then if( trim(var) == 'U' ) match = .true. if( trim(var) == 'UGRD' ) match = .true. endif if( trim(name) == 'V' ) then if( trim(var) == 'V' ) match = .true. if( trim(var) == 'VGRD' ) match = .true. endif if( trim(name) == 'T' ) then if( trim(var) == 'T' ) match = .true. if( trim(var) == 'TMPU' ) match = .true. if( trim(var) == 'TMP' ) match = .true. endif if( trim(name) == 'QV' ) then if( trim(var) == 'Q' ) match = .true. if( trim(var) == 'QV' ) match = .true. if( trim(var) == 'SPHU' ) match = .true. endif if( trim(name) == 'TH' ) then if( trim(var) == 'TH' ) match = .true. if( trim(var) == 'THETA' ) match = .true. endif if( trim(name) == 'TV' ) then if( trim(var) == 'TV' ) match = .true. endif if( trim(name) == 'THV' ) then if( trim(var) == 'THV' ) match = .true. if( trim(var) == 'THETAV' ) match = .true. endif if( trim(name) == 'PS' ) then if( trim(var) == 'PS' ) match = .true. endif if( trim(name) == 'TS' ) then if( trim(var) == 'TS' ) match = .true. endif if( trim(name) == 'O3' ) then if( trim(var) == 'O3' ) match = .true. if( trim(var) == 'OZONE' ) match = .true. endif if( trim(name) == 'DP' ) then if( trim(var) == 'DP' ) match = .true. if( trim(var) == 'DELP' ) match = .true. endif if( trim(name) == 'PHIS' ) then if( trim(var) == 'PHIS' ) match = .true. if( trim(var) == 'GZ' ) match = .true. endif return end function match end module GEOS_mkiauGridCompMod