! +-======-+ ! 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.16.2.3.2.2 2013-05-14 16:23:34 ltakacs Exp $ #define debug 0 #include "MAPL_Generic.h" !============================================================================= !BOP ! !MODULE: GEOS_mkiau -- A Module to compute the IAU forcing ! !INTERFACE: module GEOS_mkiauGridCompMod ! !USES: use ESMF use MAPL_Mod use G3_MPI_Util_Mod use GEOS_UtilsMod 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_METHOD_RUN, 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(:,:,:) :: t_glo_ana real, allocatable, dimension(:,:,:) :: q_glo_ana real, allocatable, dimension(:,:,:) :: t_glo_bkg real, allocatable, dimension(:,:,:) :: q_glo_bkg real, allocatable, dimension(:,:,:) :: r_glo_ana real, allocatable, dimension(:,:,:) :: r_glo_bkg real, allocatable, dimension(:,:,:) :: t_not real, allocatable, dimension(:,:,:) :: q_not real, allocatable, dimension(:,:,:) :: delt real, allocatable, dimension(:,:,:) :: delq real, allocatable, dimension(:,:,:) :: dp_bkg real, allocatable, dimension(:,:,:) :: rh_not real, allocatable, dimension(:,:,:) :: rh_bkg real, allocatable, dimension(:,:,:) :: rh_ana real, allocatable, dimension(:,:,:) :: ple_ana real, allocatable, dimension(:,:,:) :: pk_ana real, allocatable, dimension(:,:,:) :: pke_ana real, allocatable, dimension(:,:) :: pl, qstar, gamma real, allocatable, dimension(:,:) :: latz real, allocatable, dimension(:,:) :: g1, g2, gg 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(:,:) real, pointer :: LATS (:,:) type(ESMF_FieldBundle) :: bundle type(ESMF_Grid) :: grid type(ESMF_Time) :: currtime real :: DAMPBEG, DAMPEND, FAC real :: alf, bet, LATC real :: qave_ana, qave_bkg, qave_ana2 real :: pave_ana, pave_bkg, lambda integer :: i,j,L,n integer :: NX,NY,IMG,JMG integer :: method integer :: OFFLINE logical :: doremap logical :: FOUND logical :: ANALYZE_TS character(len=ESMF_MAXSTR) :: REPLAY_U, REPLAY_V, REPLAY_T, REPLAY_QV, REPLAY_TS, REPLAY_O3, REPLAY_RH_BKG character(len=ESMF_MAXSTR) :: REPLAY_P, REPLAY_PS, REPLAY_DP character(len=ESMF_MAXSTR) :: REPLAY_PHIS character(len=ESMF_MAXSTR) :: REPLAY_T_TYPE integer :: REPLAY_Q_CASE logical :: L_REPLAY_P, L_REPLAY_U, L_REPLAY_V, L_REPLAY_T, L_REPLAY_QV, L_REPLAY_TS, L_REPLAY_O3, L_REPLAY_RH_BKG 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, & LATS=LATS, & 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) if( ANALYZE_TS ) then REPLAY_TS = 'YES' else REPLAY_TS = 'NO' endif call MAPL_GetResource(MAPL, REPLAY_PHIS, Label="REPLAY_PHIS:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_TS, Label="REPLAY_TS:", default=trim(REPLAY_TS), RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_P , Label="REPLAY_P:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_PS, Label="REPLAY_PS:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_DP, Label="REPLAY_DP:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_U , Label="REPLAY_U:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_V , Label="REPLAY_V:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_T , Label="REPLAY_T:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_QV, Label="REPLAY_QV:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_O3, Label="REPLAY_O3:", default='YES', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_RH_BKG, Label="REPLAY_RH_BKG:", default='NO', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_T_TYPE, Label="REPLAY_T_TYPE:", default='NULL', RC=STATUS) VERIFY_(STATUS) call MAPL_GetResource(MAPL, REPLAY_Q_CASE, Label="REPLAY_Q_CASE:", default=0, RC=STATUS) VERIFY_(STATUS) REPLAY_PHIS = uppercase(REPLAY_PHIS) REPLAY_TS = uppercase(REPLAY_TS) REPLAY_P = uppercase(REPLAY_P ) REPLAY_U = uppercase(REPLAY_U ) REPLAY_V = uppercase(REPLAY_V ) REPLAY_T = uppercase(REPLAY_T ) REPLAY_QV = uppercase(REPLAY_QV) REPLAY_O3 = uppercase(REPLAY_O3) REPLAY_RH_BKG = uppercase(REPLAY_RH_BKG) REPLAY_T_TYPE = uppercase(REPLAY_T_TYPE) L_REPLAY_TS = trim(REPLAY_TS) .ne.'NO' .and. trim(REPLAY_TS) .ne.'NULL' L_REPLAY_P = trim(REPLAY_P) .ne.'NO' .and. trim(REPLAY_P) .ne.'NULL' L_REPLAY_U = trim(REPLAY_U) .ne.'NO' .and. trim(REPLAY_U) .ne.'NULL' L_REPLAY_V = trim(REPLAY_V) .ne.'NO' .and. trim(REPLAY_V) .ne.'NULL' L_REPLAY_T = trim(REPLAY_T) .ne.'NO' .and. trim(REPLAY_T) .ne.'NULL' L_REPLAY_QV = trim(REPLAY_QV) .ne.'NO' .and. trim(REPLAY_QV) .ne.'NULL' L_REPLAY_O3 = trim(REPLAY_O3) .ne.'NO' .and. trim(REPLAY_O3) .ne.'NULL' L_REPLAY_RH_BKG = trim(REPLAY_RH_BKG).ne.'NO' .and. trim(REPLAY_RH_BKG).ne.'NULL' 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( RC=STATUS) VERIFY_(STATUS) call ESMF_FieldBundleSet(RBundle, grid=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, fieldNameList=RNAMES, rc=STATUS ) VERIFY_(STATUS) if( first ) then if(MAPL_AM_I_ROOT() ) then print * print *, 'REPLAY File Variables, NQ: ', nq print *, '--------------------------' do k=1,nq print *, k,') ',trim(rnames(k)) enddo print * print *, 'REPLAY Options:' print *, '---------------' print *, 'REPLAY_TS ....... ',trim(REPLAY_TS) print *, 'REPLAY_P ........ ',trim(REPLAY_P) print *, 'REPLAY_U ........ ',trim(REPLAY_U) print *, 'REPLAY_V ........ ',trim(REPLAY_V) print *, 'REPLAY_T ........ ',trim(REPLAY_T) print *, 'REPLAY_QV ....... ',trim(REPLAY_QV) print *, 'REPLAY_O3 ....... ',trim(REPLAY_O3) print *, 'REPLAY_RH_BKG ... ',trim(REPLAY_RH_BKG) print * endif 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 allocate ( dp_bkg(IM,JM,LM) ) do L=1,lm dp_bkg(:,:,L) = ple_bkg(:,:,L)-ple_bkg(:,:,L-1) end do ! ********************************************************************** ! **** Get Pointers to ANA Replay Data **** ! ********************************************************************** allocate ( g1(IM,JM) ) allocate ( g2(IM,JM) ) allocate ( gg(IM,JM) ) allocate ( pl(IM,JM) ) allocate ( latz(IM,JM) ) allocate ( qstar(IM,JM) ) allocate ( gamma(IM,JM) ) allocate ( delt(IM,JM, LM) ) allocate ( delq(IM,JM, LM) ) 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) ) allocate ( rh_not(IM,JM, LM) ) allocate ( rh_ana(IM,JM, LM) ) allocate ( rh_bkg(IM,JM, LM) ) allocate ( t_not(IM,JM, LM) ) allocate ( q_not(IM,JM, LM) ) doremap = trim(cremap).eq.'YES' ! Surface Temperature !-------------------- if( L_REPLAY_TS ) then FOUND = .false. do k=1,nq if( match('ts',REPLAY_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',REPLAY_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( L_REPLAY_P ) then FOUND = .false. do k=1,nq if( match('ps',REPLAY_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( L_REPLAY_P ) then FOUND = .false. do k=1,nq if( match('dp',REPLAY_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( L_REPLAY_U ) then FOUND = .false. do k=1,nq if( match('u',REPLAY_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( L_REPLAY_V ) then FOUND = .false. do k=1,nq if( match('v',REPLAY_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( L_REPLAY_QV ) then FOUND = .false. do k=1,nq if( match('qv',REPLAY_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( L_REPLAY_O3 ) then FOUND = .false. do k=1,nq if( match('o3',REPLAY_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( L_REPLAY_T ) then if( trim(REPLAY_T).ne.'YES' .and. trim(REPLAY_T_TYPE).eq.'NULL' ) then if(MAPL_AM_I_ROOT()) then print * print *, 'You must specify REPLAY_T_TYPE when setting REPLAY_T name.' print *, 'REPLAY_T_TYPE OPTIONS: T, TV, TH, THV' print * endif RETURN_(ESMF_FAILURE) endif FOUND = .false. do k=1,nq if( match('t',REPLAY_T,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then if( trim(REPLAY_T_TYPE).eq.'NULL' .or. trim(REPLAY_T_TYPE).eq.'T' ) then t_ana = TEMP3D FOUND = .true. exit endif endif endif if( match('tv',REPLAY_T,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then if( trim(REPLAY_T_TYPE).eq.'NULL' .or. trim(REPLAY_T_TYPE).eq.'TV' ) then t_ana = TEMP3D/(1.0+eps*q_ana) FOUND = .true. exit endif endif endif if( match('th',REPLAY_T,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then if( trim(REPLAY_T_TYPE).eq.'NULL' .or. trim(REPLAY_T_TYPE).eq.'TH' ) then t_ana = TEMP3D*pk_ana FOUND = .true. exit endif endif endif if( match('thv',REPLAY_T,rnames(k)) ) then call ESMFL_BundleGetPointertoData(RBundle,trim(rnames(k)),TEMP3D, RC=STATUS) if(STATUS==ESMF_SUCCESS) then if( trim(REPLAY_T_TYPE).eq.'NULL' .or. trim(REPLAY_T_TYPE).eq.'THV' ) then t_ana = TEMP3D*pk_ana/(1.0+eps*q_ana) FOUND = .true. exit endif 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 ! Construct RH_BKG ! ---------------- do L=1,lm pl(:,:) = ( ple_bkg(:,:,L)+ple_bkg(:,:,L-1) )*0.5 qstar(:,:) = GEOS_QSAT( t_bkg(:,:,L),pl(:,:),PASCALS=.true. ) rh_bkg(:,:,L) = q_bkg(:,:,L)/qstar(:,:) enddo ! 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 ! Construct RH_ANA ! ---------------- do L=1,lm pl(:,:) = ( ple_ana(:,:,L)+ple_ana(:,:,L-1) )*0.5 qstar(:,:) = GEOS_QSAT( t_ana(:,:,L),pl(:,:),PASCALS=.true. ) rh_ana(:,:,L) = q_ana(:,:,L)/qstar(:,:) enddo ! Get Global Dimensions and LAYOUT ! -------------------------------- 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) ! Adjust T_ANA and Q_ANA to Replay RH_BKG ! --------------------------------------- if( L_REPLAY_RH_BKG ) then #if debug call create_dynamics_lattice ( lattice,mp_comm_world,IMG,JMG,lm,NX,NY) allocate( t_glo_ana(img,jmg,lm) ) allocate( q_glo_ana(img,jmg,lm) ) allocate( r_glo_ana(img,jmg,lm) ) allocate( t_glo_bkg(img,jmg,lm) ) allocate( q_glo_bkg(img,jmg,lm) ) allocate( r_glo_bkg(img,jmg,lm) ) call G3_GATHER ( t_glo_ana, t_ana,lattice ) call G3_GATHER ( q_glo_ana, q_ana,lattice ) call G3_GATHER ( r_glo_ana,rh_ana,lattice ) call G3_GATHER ( t_glo_bkg, t_bkg,lattice ) call G3_GATHER ( q_glo_bkg, q_bkg,lattice ) call G3_GATHER ( r_glo_bkg,rh_bkg,lattice ) if(MAPL_AM_I_ROOT() ) then print *, "writing initial data ..." do L=lm,1,-1 write(68) t_glo_ana(:,:,L) enddo do L=lm,1,-1 write(68) q_glo_ana(:,:,L)*1000 enddo do L=lm,1,-1 write(68) r_glo_ana(:,:,L)*100 enddo do L=lm,1,-1 write(68) t_glo_bkg(:,:,L) enddo do L=lm,1,-1 write(68) q_glo_bkg(:,:,L)*1000 enddo do L=lm,1,-1 write(68) r_glo_bkg(:,:,L)*100 enddo endif #endif if( trim(REPLAY_RH_BKG).eq.'YES' ) then ! Minimize ( Cp*DelT )**2 and (L*DelQ)**2 alf = MAPL_CP**2 bet = MAPL_ALHL**2 endif if( trim(REPLAY_RH_BKG).eq.'RHE' ) then ! Minimize ( DelE )**2 where E = Cp*T + L*Q alf = MAPL_CP bet = MAPL_ALHL endif FAC = 0.622*MAPL_ALHL/MAPL_RGAS t_not = t_ana q_not = q_ana do L=1,lm do j=1,jm do i=1,im pl(i,j) = ( ple_ana(i,j,L)+ple_ana(i,j,L-1) )*0.5 do n=1,10 qstar(i,j) = GEOS_QSAT( t_not(i,j,L),pl(i,j),PASCALS=.true. ) rh_not(i,j,L) = q_ana(i,j,L)/qstar(i,j) gamma(i,j) = FAC/t_not(i,j,L)**2 ! Minimize ( Cp*DelT )**2 and (L*DelQ)**2 ! --------------------------------------- if( trim(REPLAY_RH_BKG).eq.'YES' ) then delt(i,j,L) = - ( alf * ( t_not(i,j,L)-t_ana(i,j,L) ) & + bet*rh_bkg(i,j,L)*gamma(i,j)*qstar(i,j)**2 *( rh_bkg(i,j,L)-rh_not(i,j,L) ) ) & / ( alf + bet*rh_bkg(i,j,L)*gamma(i,j)*qstar(i,j)**2 *( (rh_bkg(i,j,L)-rh_not(i,j,L) ) * & (gamma(i,j)-2/t_not(i,j,L))+ rh_bkg(i,j,L)*gamma(i,j) ) ) endif ! Minimize ( DelE )**2 where E = Cp*T + L*Q ! ----------------------------------------- if( trim(REPLAY_RH_BKG).eq.'RHE' ) then delt(i,j,L) = - ( alf*( t_not(i,j,L)-t_ana(i,j,L) ) + bet*qstar(i,j)*( rh_bkg(i,j,L)-rh_not(i,j,L) ) ) & / ( alf + bet*rh_bkg(i,j,L)*gamma(i,j)*qstar(i,j) ) endif t_not(i,j,L) = t_not(i,j,L) + delt(i,j,L) enddo qstar(i,j) = GEOS_QSAT( t_not(i,j,L),pl(i,j),PASCALS=.true. ) q_not(i,j,L) = qstar(i,j)*rh_bkg(i,j,L) enddo enddo enddo ! Blend Resulting T_NOT and Q_NOT with Initial T_ANA and Q_ANA Values ! ------------------------------------------------------------------- LATC = 20.0 LATZ = abs(LATS*180.0/MAPL_PI)-LATC G1 = 1.0-0.7*EXP(-abs(LATZ)/5.0) G2 = 0.3*EXP(-abs(LATZ)/5.0) WHERE ( LATZ.ge.0.0 ) GG = G1 ELSEWHERE GG = G2 ENDWHERE do L=1,lm do j=1,jm do i=1,im t_ana(i,j,L) = gg(i,j)*t_not(i,j,L) + (1-gg(i,j))*t_ana(i,j,L) q_ana(i,j,L) = gg(i,j)*q_not(i,j,L) + (1-gg(i,j))*q_ana(i,j,L) pl(i,j) = ( ple_ana(i,j,L)+ple_ana(i,j,L-1) )*0.5 qstar(i,j) = GEOS_QSAT( t_ana(i,j,L),pl(i,j),PASCALS=.true. ) rh_ana(i,j,L) = q_ana(i,j,L)/qstar(i,j) enddo enddo enddo #if debug call G3_GATHER ( t_glo_ana, t_ana,lattice ) call G3_GATHER ( q_glo_ana, q_ana,lattice ) call G3_GATHER ( r_glo_ana,rh_ana,lattice ) if(MAPL_AM_I_ROOT() ) then print *, "writing final data ..." do L=lm,1,-1 write(68) t_glo_ana(:,:,L) enddo do L=lm,1,-1 write(68) q_glo_ana(:,:,L)*1000 enddo do L=lm,1,-1 write(68) r_glo_ana(:,:,L)*100 enddo do L=lm,1,-1 write(68) t_glo_bkg(:,:,L) enddo do L=lm,1,-1 write(68) q_glo_bkg(:,:,L)*1000 enddo do L=lm,1,-1 write(68) r_glo_bkg(:,:,L)*100 enddo endif deallocate( t_glo_ana ) deallocate( q_glo_ana ) deallocate( r_glo_ana ) deallocate( t_glo_bkg ) deallocate( q_glo_bkg ) deallocate( r_glo_bkg ) call destroy_dynamics_lattice (lattice) #endif endif ! End REPLAY_RH_BKG Test #if 0 ! Blend Q_BKG and Q_ANA Values ! ---------------------------- LATC = 20.0 LATZ = abs(LATS*180.0/MAPL_PI)-LATC G1 = 1.0-0.7*EXP(-abs(LATZ)/5.0) G2 = 0.3*EXP(-abs(LATZ)/5.0) WHERE ( LATZ.ge.0.0 ) GG = G1 ELSEWHERE GG = G2 ENDWHERE do L=1,lm do j=1,jm do i=1,im q_ana(i,j,L) = gg(i,j)*q_bkg(i,j,L) + (1-gg(i,j))*q_ana(i,j,L) enddo enddo enddo #endif ! ********************************************************************** ! **** Blend ANA and BKG Variables between DAMPBEG and DAMPEND **** ! ********************************************************************** call create_dynamics_lattice ( lattice,mp_comm_world,IMG,JMG,lm,NX,NY) 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) ! ------------------------------------------------------------------------ if( JMG.ne.6*IMG ) then method = 1 call windfix ( u_ana,v_ana,ple_ana, & u_bkg,v_bkg,ple_bkg,im,jm,lm,lattice,'A',method ) endif endif ! Perform Global Constraint on Vertically Integrated Moisture Increment ! --------------------------------------------------------------------- if( REPLAY_Q_CASE .ne. 0 ) then call gmean ( q_ana,qave_ana ,dp_ana,im,jm,lm,1,lattice ) call gmean ( q_bkg,qave_bkg ,dp_bkg,im,jm,lm,1,lattice ) if( REPLAY_Q_CASE .eq. 1 ) q_not = q_ana if( REPLAY_Q_CASE .eq. 2 ) q_not = q_ana-q_bkg call gmean ( q_not,qave_ana2,dp_ana,im,jm,lm,2,lattice ) q_not = 1.0 call gmean ( q_not,pave_ana ,dp_ana,im,jm,lm,1,lattice ) call gmean ( q_not,pave_bkg ,dp_bkg,im,jm,lm,1,lattice ) if(MAPL_AM_I_ROOT() ) then print * print *, 'P Before: ',pave_ana/100,pave_bkg/100,(pave_ana-pave_bkg)/100 print *, 'Q Before: ',qave_ana/100,qave_bkg/100,(qave_ana-qave_bkg)/100 endif ! lambda = ( qave_ana-qave_bkg-pave_ana+pave_bkg )/qave_ana2 lambda = ( qave_ana-qave_bkg )/qave_ana2 q_ana = q_ana - (q_ana-q_bkg)**2 * lambda call gmean ( q_ana,qave_ana ,dp_ana,im,jm,lm,1,lattice ) if(MAPL_AM_I_ROOT() ) then print *, 'Q After: ',qave_ana/100,qave_bkg/100,(qave_ana-qave_bkg)/100 print *, ' lambda: ',lambda endif endif call destroy_dynamics_lattice (lattice) ! ********************************************************************** ! **** 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( L_REPLAY_U ) then du = u_ana-u_bkg else du = 0.0 endif if( L_REPLAY_V ) then dv = v_ana-v_bkg else dv = 0.0 endif if( L_REPLAY_T ) then dt = t_ana-t_bkg else dt = 0.0 endif if( L_REPLAY_QV ) then dq = q_ana-q_bkg else dq = 0.0 endif if( L_REPLAY_O3 ) then do3 = o3_ana-o3_bkg else do3 = 0.0 endif if( L_REPLAY_P ) then dple = ple_ana-ple_bkg else dple = 0.0 endif if( L_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_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 ( dp_bkg ) deallocate ( pk_ana ) deallocate ( pke_ana ) deallocate ( thv_ana ) deallocate ( rh_not ) deallocate ( rh_ana ) deallocate ( rh_bkg ) deallocate ( t_not ) deallocate ( q_not ) deallocate ( latz ) deallocate ( qstar ) deallocate ( gamma ) deallocate ( delt ) deallocate ( delq ) deallocate ( pl,gg ) deallocate ( g1,g2 ) 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