#define VERIFY_(rc) If (rc /= ESMF_SUCCESS) Call abort_core(__LINE__,rc) #define RETURN_(status) If (Present(rc)) rc=status; return #if defined(USE_FVCORE) !#define NO_FORCING !--------------------------------------------------------------------------------------------- ! The following module provides an interface for modelE to use the ESMF-wrapped version of ! the FV dynamical core from GEOS-5. Many elements closely resemble an ESMF coupler component, ! which is to be expected given the intended role for this layer. !--------------------------------------------------------------------------------------------- module FV_INTERFACE_MOD 8,10 use esmf_mod !!$ use GEOS_Mod, Only: ESMFL_StateGetPointerToData use ESMFL_MOD, Only: ESMFL_StateGetPointerToData implicit none private ! except for public :: Write_Profile public :: FV_CORE ! Derived type to encapsulate FV + modelE data public :: Init_app_clock ! modelE does not use ESMF clocks, but ESMF requires one public :: Initialize ! Uses modelE data to initialize the FV gridded component public :: Finalize ! Uses modelE data to finalize the FV gridded component public :: Compute_Tendencies ! Temporarily a separate method, but will move within Run() soon. public :: Run ! Execute the FV method to integrate the core forward in time public :: Checkpoint ! Unimplemented public :: edgepressure_giss public :: DryTemp_GISS ! This data structure is o convenient entity for storing persistent data between ! calls to this module. In addition to Type FV_CORE PRIVATE type (esmf_gridcomp) :: gc ! This is the handle for the fv dynamical core type (esmf_grid) :: grid ! Although modelE is not an ESMF component, it does have an ESMF_Grid type (esmf_vm) :: vm ! Should be eliminated ... Only used to obtain NPES for config file. ! Import and Export states for FV dycore type(esmf_state) :: import ! Allocated within FV component type(esmf_state) :: export ! Allocated within FV component ! The following pointers can be re-extracted from the import state at each iteration, ! but it is convenient to have a simpler means of access. real*4, pointer, dimension(:,:,:) :: dudt, dvdt, dtdt, dpedt ! Tendencies real*4, pointer, dimension(:,:,:) :: Q ! Humidity real*4, pointer, dimension(:,:,:) :: Qtr ! other tracers real*4, pointer, dimension(:,:) :: phis ! modelE does not work directly with tendencies. Instead, tendencies are derived ! by differencing before and after physics. Therefore, the final dynamical state ! must be preserved for the following ! of modelE fields real*8, pointer, dimension(:,:,:) :: U_old, V_old, dPT_old, PE_old, dT_old END Type FV_CORE ! private data of convenience Integer :: rc ! return code from ESMF Interface Initialize 1,2 module procedure initialize_fv End Interface Interface Run 1 module procedure run_fv End Interface Interface reverse Module Procedure reverse_3d_r8 Module Procedure reverse_3d_r4 End Interface Interface ConvertPotTemp_GISS2FV 1 Module Procedure CnvPotTemp_GISS2FV_r8 Module Procedure CnvPotTemp_GISS2FV_r4 End Interface Interface ConvertPotTemp_FV2GISS 1 Module Procedure CnvPotTemp_FV2GISS_r4 End Interface Interface ConvertPressure_GISS2FV 2 Module Procedure ConvertPressure_GISS2FV_r4 Module Procedure ConvertPressure_GISS2FV_r8 End Interface ! The following parameters address the fact that FV and modelE use ! different units and different reference pressures for potential temperature. ! Superficially there is redundancy between the two sets, but in some sense ! this is merely coincidental. Real*8, parameter :: REF_PRESSURE_GISS = 100 ! 1 mb = 100 pa Real*8, parameter :: REF_PRESSURE_FV = 1 ! pa Real*8, parameter :: REF_RATIO = REF_PRESSURE_FV / REF_PRESSURE_GISS Real*8, parameter :: PRESSURE_UNIT_GISS = 100 ! 1 mb Real*8, parameter :: PRESSURE_UNIT_FV = 1 ! 1 pa Real*8, parameter :: PRESSURE_UNIT_RATIO = PRESSURE_UNIT_GISS/PRESSURE_UNIT_FV character(len=*), parameter :: FVCORE_INTERNAL_RESTART = 'fv_internal_restart.dat' character(len=*), parameter :: FVCORE_IMPORT_RESTART = 'fv_import_restart.dat' character(len=*), parameter :: FVCORE_LAYOUT = 'fv_layout.rc' character(len=*), parameter :: TENDENCIES_FILE = 'tendencies_checkpoint' integer, parameter :: initial_start = 8 integer, parameter :: extend_run = 9 contains subroutine Initialize_fv(fv, istart, vm, grid, clock, config_file) 1,14 use resolution, only: IM, LM use fvdycore_gridcompmod, only: SetServices use GEOS_BaseMod, only: GEOS_FieldCreate use DOMAIN_DECOMP, only : modelE_grid => grid, get type (fv_core), intent(inout) :: fv integer, intent(in) :: istart type (esmf_vm), intent(in) :: vm type (esmf_grid), intent(in) :: grid type (esmf_clock), intent(in) :: clock character(len=*), intent(in) :: config_file ! filename for resource file type (esmf_config) :: cf Integer :: L Integer :: j_0, j_1 type (ESMF_Bundle) :: bundle type (ESMF_Field) :: Qfield type (ESMF_Array) :: Qarray type (ESMF_FieldDataMap) :: datamap fv % vm =vm fv % grid=grid ! Load configuration information from resource file ! ------------------------------------------------ cf = load_configuration(config_file) ! Create the dynamics component, using same layout as application ! ------------------------------------------------------------------ fv % gc = ESMF_GridCompCreate ( vm=vm, name='FV dynamics', & & grid=grid, gridcomptype=ESMF_ATM, & & config=cf, rc=rc) VERIFY_(rc) ! Create couplings fv % import = ESMF_StateCreate ( 'fv dycore imports', ESMF_STATE_IMPORT, rc=rc ) VERIFY_(rc) fv % export = ESMF_StateCreate ( 'fv dycore exports', ESMF_STATE_EXPORT, rc=rc ) VERIFY_(rc) ! Register services for components ! -------------------------------- ! print*,'calling set services' call ESMF_GridCompSetServices ( fv % gc, SetServices, rc ) VERIFY_(rc) ! The FV components requires its own restart file for managing ! its internal state. We check to see if the file already exists, and if not ! create one based upon modelE's internal state. call create_restart_file(fv, istart, cf, clock) ! print*,'calling fv initialize' call ESMF_GridCompInitialize ( fv % gc, importState=fv % import, exportState=fv % export, clock=clock, & & phase=ESMF_SINGLEPHASE, rc=rc ) VERIFY_(rc) ! Initialize component (and import/export states) ! ---------------------------------------------- Call allocate_tendency_storage(fv, istart) call allocateFvExport3D ( fv % export,'U' ) call allocateFvExport3D ( fv % export,'V' ) call allocateFvExport3D ( fv % export,'TH' ) call allocateFvExport3D ( fv % export,'PLE' ) call allocateFvExport3D ( fv % export,'Q' ) call allocateFvExport3D ( fv % export,'MFX' ) call allocateFvExport3D ( fv % export,'MFY' ) call allocateFvExport3D ( fv % export,'MFZ' ) ! print*,'called fv initialize' ! Specific Humidity - need to reserve space in FV import state call ESMFL_StateGetPointerToData ( fv % export,fv % q,'Q',rc=rc) VERIFY_(rc) Call get(modelE_grid, J_STRT=J_0, J_STOP=J_1) !!$ Allocate(fv % Qtr(IM,J_0:J_1,LM)) Qarray = ESMF_ArrayCreate(fv % q, ESMF_DATA_REF, RC=rc) VERIFY_(rc) call ESMF_FieldDataMapSetDefault(datamap, dataRank=3, rc=rc) Qfield = ESMF_FieldCreate( grid, Qarray, horzRelloc=ESMF_CELL_CENTER, & datamap=datamap, name="Q", rc=rc) VERIFY_(rc) ! First obtain a reference field from the state (avoids tedious details) call ESMF_StateGetBundle (fv % import, 'QTR' , bundle, rc=rc) VERIFY_(rc) Call ESMF_BundleAddField(bundle, Qfield, rc=rc) VERIFY_(rc) contains subroutine allocateFvExport3D ( state, name ) 8 type(ESMF_State), intent(INOUT) :: state character(len=*), intent(IN ) :: name real, pointer :: ptr(:,:,:) logical :: alloc integer :: status alloc = .true. call ESMFL_StateGetPointerToData ( state, ptr , name , alloc , status ) VERIFY_(status) end subroutine allocateFvExport3D end subroutine Initialize_fv function EdgePressure_GISS() Result(PE) 5,5 USE RESOLUTION, only: IM, LM, LS1 Use MODEL_COM, only : SIG, SIGE, Ptop, PSFMPT, P use DOMAIN_DECOMP, only: grid, get USE CONSTANT, only: KAPA REAL*8 :: PE(IM,grid % J_STRT:grid % J_STOP,LM+1) INTEGER :: L, j_0, j_1 call get(grid, J_STRT=J_0, J_STOP=J_1) Do L = 1, LM+1 If (L < LS1) THEN PE(:,:,L) = SIGE(L)*P(:,J_0:J_1) + Ptop Else PE(:,:,L) = SIGE(L)*PSFMPT + Ptop End IF End Do end function EdgePressure_GISS ! Compute Delta-pressure for GISS model function DeltPressure_GISS() Result(dP) 1,5 USE RESOLUTION, only: IM, LM, LS1 Use MODEL_COM, only: T USE DOMAIN_DECOMP, only: grid, GET REAL*8 :: dP(IM,grid % J_STRT:grid % J_STOP,LM) REAL*8 :: PE(IM,grid % J_STRT:grid % J_STOP,LM+1) INTEGER :: J_0,J_1,k Call Get(grid, J_STRT=J_0, J_STOP=J_1) PE = EdgePressure_GISS() do k = 1, LM dP(:,:,k) = (PE(:,:,k)-PE(:,:,k+1)) end do end function DeltPressure_GISS ! Convert Potential Temperature into (dry) Temperature function DeltPressure_DryTemp_GISS() Result(dPT) 2,6 USE RESOLUTION, only: IM, LM, LS1 Use MODEL_COM, only: T USE DOMAIN_DECOMP, only: grid, GET REAL*8 :: dPT(IM,grid % J_STRT:grid % J_STOP,LM) REAL*8 :: PE(IM,grid % J_STRT:grid % J_STOP,LM+1) REAL*8 :: T_dry(IM,grid % J_STRT:grid % J_STOP,LM) INTEGER :: J_0,J_1,k Call Get(grid, J_STRT=J_0, J_STOP=J_1) T_dry = DryTemp_GISS() PE = EdgePressure_GISS() do k = 1, LM dPT(:,:,k) = (PE(:,:,k)-PE(:,:,k+1)) * T_dry(:,:,k) end do end function DeltPressure_DryTemp_GISS Subroutine allocate_tendency_storage(fv, istart) 1,16 Use DOMAIN_DECOMP, only: GRID, GET, AM_I_ROOT USE RESOLUTION, only: IM, LM, LS1 USE MODEL_COM, only: U, V, T use FILEMANAGER Implicit None Type (FV_Core) :: fv integer, intent(in) :: istart Integer :: J_0H, J_1H, J_0, J_1 integer :: iunit Call Get(grid, J_strt_halo=J_0H, J_stop_halo=J_1H, & & J_STRT=J_0, J_STOP=J_1) ! 1) Link/copy modelE data to import state call ESMFL_StateGetPointerToData ( fv % import,fv % dudt,'DUDT',rc=rc) VERIFY_(rc) call ESMFL_StateGetPointerToData ( fv % import,fv % dvdt,'DVDT',rc=rc) VERIFY_(rc) call ESMFL_StateGetPointerToData ( fv % import,fv % dtdt,'DTDT',rc=rc) VERIFY_(rc) call ESMFL_StateGetPointerToData ( fv % import,fv % dpedt,'DPEDT',rc=rc) VERIFY_(rc) call ESMFL_StateGetPointerToData ( fv % import,fv % phis, 'PHIS', rc=rc) VERIFY_(rc) ! 2) Allocate space for storing old values from which to compute tendencies Allocate(fv % U_old(IM,J_0:J_1,LM), & & fv % V_old(IM,J_0:J_1,LM), & & fv % dPT_old(IM,J_0:J_1,LM), & & fv % dT_old(IM,J_0:J_1,LM), & & fv % PE_old(IM,J_0:J_1,LM+1)) select case (istart) case (:initial_start) ! Do a cold start. Set Old = Current. fv % dudt=0 fv % dvdt=0 fv % dtdt=0 fv % dpedt=0 fv % U_old = U(:,J_0:J_1,:) fv % V_old = V(:,J_0:J_1,:) fv % dPT_old = DeltPressure_DryTemp_GISS() fv % dT_old = DryTemp_GISS() fv % PE_old = EdgePressure_GISS() case (extend_run:) ! input couplings are in a file somewhere and already read in if (AM_I_ROOT()) call openunit(TENDENCIES_FILE, iunit, qbin=.true.,qold=.true.) call readArr(iunit, fv % U_old) call readArr(iunit, fv % V_old) call readArr(iunit, fv % dPT_old) call readArr(iunit, fv % PE_old) call readArr(iunit, fv % dT_old) if (AM_I_ROOT()) call closeunit(iunit) call compute_tendencies(fv) !! case default !! call stop_model('ISTART option not supported',istart) end select contains subroutine readArr(iunit, arr) 5,3 use resolution, only: IM, JM use DOMAIN_DECOMP, only: grid, unpack_data, get integer, intent(in) :: iunit real*8, intent(out) :: arr(:,:,:) real*8, allocatable :: padArr(:,:,:) real*8, allocatable :: globalArr(:,:,:) allocate(globalArr(IM,JM,size(arr,3))) allocate(padArr(1:IM,j_0h:j_1h, size(arr,3))) if (AM_I_ROOT()) read(iunit) globalArr call unpack_data(grid, globalArr, padArr) arr(:,:,:) = padArr(:,j_0:j_1,:) deallocate(padArr) deallocate(globalArr) end subroutine readArr End Subroutine allocate_tendency_storage subroutine finalize(fv, clock, fv_fname, fv_dfname) 1,3 USE DOMAIN_DECOMP, only: am_I_root Type (FV_Core) :: fv type (esmf_clock), intent(in) :: clock integer :: rc character(len=*), intent(in) :: fv_fname, fv_dfname call saveTendencies(fv, fv_dfname) call ESMF_GridCompFinalize ( fv % gc, fv % import, fv % export, clock, rc=rc ) if (AM_I_ROOT()) then call system('mv ' // FVCORE_INTERNAL_RESTART // ' ' // trim(fv_fname) ) call system('rm ' // TENDENCIES_FILE ) end if Deallocate(fv % U_old, fv % V_old, fv % dPT_old, fv % dT_old, fv % PE_old) end subroutine finalize subroutine clearTendencies(fv) 2 Implicit None Type (FV_CORE) :: fv fv % dudt = 0 fv % dvdt = 0 fv % dtdt = 0 fv % dpedt = 0 end subroutine clearTendencies ! Compute tendencies ! ------------------ subroutine compute_tendencies(fv) 2,16 USE RESOLUTION, only: IM, LM USE MODEL_COM, Only: U, V, T USE DOMAIN_DECOMP, only: grid, get USE DYNAMICS, only: DUT, DVT USE CONSTANT, only: KAPA Implicit None Type (FV_CORE) :: fv Integer :: J_0, J_1 Integer :: J_0h, J_1h Call Get(grid, j_strt=j_0, j_stop=j_1) ! U, V DUT(:,J_0:J_1,:) = Tendency(U(:,J_0:J_1,:), fv % U_old(:,J_0:J_1,:)) DVT(:,J_0:J_1,:) = Tendency(V(:,J_0:J_1,:), fv % V_old(:,J_0:J_1,:)) call ConvertUV_GISS2FV(DUT(:,J_0:J_1,:), DVT(:,J_0:J_1,:), fv % dudt, fv % dvdt) ! delta pressure weighted Temperature fv % dtdt = reverse(DeltPressure_GISS() * Tendency(DryTemp_GISS(), fv % dT_old)) * & & (PRESSURE_UNIT_RATIO) ! Edge Pressure Call ConvertPressure_GISS2FV( Tendency(EdgePressure_GISS(), fv % PE_old), fv % dpedt) #ifdef NO_FORCING call clearTendencies() #endif end subroutine compute_tendencies !----------- ! d/dt A ! Note that FV needs tendency with the order of vertical levels reversed. !----------- Function Tendency(A, A_old) Result(tend) 4,1 USE MODEL_COM, Only: DTsrc REAL*8, INTENT(IN) :: A(:,:,:) REAL*8, INTENT(IN) :: A_old(:,:,:) REAL*8 :: tend(size(a,1),size(a,2),size(a,3)) tend = (A - A_old)/DTsrc End Function Tendency subroutine run_fv(fv, clock) 1,19 USE DOMAIN_DECOMP, only: grid, halo_update, get, grid, NORTH USE MODEL_COM, Only : U, V, T, P, IM, JM, LM, ZATMO USE MODEL_COM, only : NIdyn, DT, DTSRC USE SOMTQ_COM, only: QMOM, TMOM, MZ USE ATMDYN, only: CALC_AMP, CALC_PIJL, AFLUX, COMPUTE_MASS_FLUX_DIAGS USE DYNAMICS, only: MA, PHI, GZ USE DYNAMICS, ONLY: PU, PV, CONV USE DYNAMICS, ONLY: SD, AFLUX, PUA, PVA, SDA Type (FV_CORE) :: fv Type (ESMF_Clock) :: clock type (ESMF_TimeInterval) :: timeInterval REAL*8, DIMENSION(IM,grid % J_STRT_HALO:grid % J_STOP_HALO,LM) :: PIJL integer :: L,istep, NS, NIdyn_fv integer :: rc !@sum CALC_AMP Calc. AMP: kg air*grav/100, incl. const. pressure strat call calc_amp(P, MA) CALL CALC_PIJL(LM,P,PIJL) Call Copy_modelE_to_FV_import(fv) call clear_accumulated_mass_fluxes() ! Run dycore NIdyn_fv = DTsrc / (DT) do istep = 1, NIdyn_fv call ESMF_GridCompRun ( fv % gc, fv % import, fv % export, clock, 91, rc=rc ) call clearTendencies(fv) call ESMF_GridCompRun ( fv % gc, fv % import, fv % export, clock, rc=rc ) call ESMF_TimeIntervalSet(timeInterval, s = nint(DT), rc=rc) call ESMF_ClockAdvance(clock, timeInterval, rc=rc) call accumulate_mass_fluxes(fv) call Copy_FV_export_to_modelE(fv) ! inside loop to accumulate PUA,PVA,SDA call reset_tmom #if defined(USE_FV_Q) call reset_qmom #endif phi = compute_phi(P, T, TMOM(MZ,:,:,:), ZATMO) call compute_mass_flux_diags(phi, pu, pv, dt) end do gz = phi end subroutine run_fv subroutine checkpoint(fv, clock, fv_fname, fv_dfname) 2,2 use GEOS_mod, only: GEOS_RecordPhase use DOMAIN_DECOMP, only: AM_I_ROOT Type (FV_Core), intent(inout) :: fv type (esmf_clock), intent(in) :: clock integer :: rc character(len=*), parameter :: SUFFIX_TEMPLATE = '.YYYYMMDD_HHMMz.bin' character(len=len(SUFFIX_TEMPLATE)) :: suffix Type (ESMF_Time) :: currentTime integer :: year, month, day, hour, minute, second character(len=*), intent(in) :: fv_fname, fv_dfname call ESMF_GridCompFinalize( fv % gc, fv % import, fv % export, clock, GEOS_RecordPhase, rc=rc) ! Now move the file into a more useful name if (AM_I_ROOT()) then ! 1) FV names restarts as: fvcore_internal_checkoint.YYYYMMDD_HHMMz.bin call ESMF_ClockGet(clock, currTime=currentTime, rc=rc) call ESMF_TimeGet(currentTime, YY=year, MM= month, & DD=day, H=hour, M=minute, & S=second, rc=rc) write(suffix,'(".",i4.4,i2.2,i2.2,"_",i2.2,i2.2,"z.bin")'), & & year, month, day, hour, minute call system('mv ' // FVCORE_INTERNAL_RESTART // suffix // ' ' // trim(fv_fname) ) end if call saveTendencies(fv, fv_dfname) end subroutine checkpoint subroutine saveTendencies(fv, fv_dfname) 2,9 use DOMAIN_DECOMP, only: AM_I_ROOT use FILEMANAGER type (fv_core), intent(inout) :: fv integer :: iunit character(len=*), intent(in) :: fv_dfname if (AM_I_ROOT()) call openunit(trim(fv_dfname) , iunit, qbin=.true.) call saveArr(iunit, fv % U_old) call saveArr(iunit, fv % V_old) call saveArr(iunit, fv % dPT_old) call saveArr(iunit, fv % PE_old) call saveArr(iunit, fv % dT_old) if (AM_I_ROOT()) call closeunit(iunit) contains subroutine saveArr(iunit, arr) 5,4 use resolution, only: IM, JM use DOMAIN_DECOMP, only: grid, pack_data, get integer, intent(in) :: iunit real*8, intent(in) :: arr(:,:,:) real*8, allocatable :: padArr(:,:,:) real*8, allocatable :: globalArr(:,:,:) integer :: j_0, j_1 integer :: j_0h, j_1h Call Get(grid, j_strt=j_0, j_stop=j_1, & & j_strt_halo = j_0h, j_stop_halo = j_1h) allocate(globalArr(IM,JM,size(arr,3))) allocate(padArr(1:IM,j_0h:j_1h, size(arr,3))) padArr(:,j_0:j_1,:) = arr(:,:,:) call pack_data(grid, padArr, globalArr) if (AM_I_ROOT()) write(iunit) globalArr deallocate(padArr) deallocate(globalArr) end subroutine saveArr end subroutine saveTendencies !---------------------------- ! Internal routines !---------------------------- !---------------------------- function init_app_clock(start_time, end_time, interval) Result(clock) 1 integer :: start_time(6) integer :: end_time(6) integer :: interval type (esmf_clock) :: clock type (esmf_time) :: startTime type (esmf_time) :: stopTime type (esmf_timeinterval) :: timeStep type (esmf_calendar) :: gregorianCalendar ! initialize calendar to be Gregorian type gregorianCalendar = esmf_calendarcreate("GregorianCalendar", ESMF_CAL_GREGORIAN, rc) VERIFY_(rc) ! initialize start time write(*,*)'Time Set Start: ',START_TIME call esmf_timeset(startTime, YY=START_TIME(1), MM= START_TIME(2), & & DD=START_TIME(3), H=START_TIME(4), & & M=START_TIME(5), S=START_TIME(6), & & calendar=gregorianCalendar, rc=rc) VERIFY_(rc) ! initialize stop time write(*,*)'Time Set End: ',END_TIME call esmf_timeset(stopTime, YY=END_TIME(1), MM= END_TIME(2), & & DD=END_TIME(3), H=END_TIME(4), & & M=END_TIME(5), S=END_TIME(6), & & calendar=gregorianCalendar, rc=rc) VERIFY_(rc) ! initialize time interval call esmf_timeintervalset(timeStep, & & S=INT(interval), rc=rc) VERIFY_(rc) ! initialize the clock with the above values clock = esmf_clockcreate("ApplClock", timeStep, startTime, stopTime, rc=rc) VERIFY_(rc) end function init_app_clock !---------------------------- !---------------------------- function load_configuration(config_file) result( config ) 1,5 use ESMF_mod use FILEMANAGER Use MODEL_COM, only: DT character(len=*), parameter :: Iam="FV_INTERFACE::loadconfiguration" character(len=*), intent(in) :: config_file type (esmf_config) :: config integer :: iunit type (ESMF_VM) :: vm config = esmf_configcreate(rc=rc) VERIFY_(rc) call openunit(config_file, iunit, qbin=.false., qold=.false.) write(iunit,*)'FVCORE_INTERNAL_CHECKPOINT_FILE: ', FVCORE_INTERNAL_RESTART write(iunit,*)'FVCORE_INTERNAL_RESTART_FILE: ', FVCORE_INTERNAL_RESTART !!$ write(iunit,*)'FVCORE_IMPORT_CHECKPOINT_FILE: ', TENDENCIES_FILE !!$ write(iunit,*)'FVCORE_IMPORT_RESTART_FILE: ', TENDENCIES_FILE write(iunit,*)'FVCORE_LAYOUT_FILE: ', FVCORE_LAYOUT write(iunit,*)'RUN_DT: ', DT call closeUnit(iunit) Call ESMF_VMGetGlobal(vm, rc) call esmF_VMbarrier(vm, rc) call esmf_configloadfile(config, config_file, rc=rc) VERIFY_(rc) end function load_configuration !---------------------------- Subroutine Create_Restart_File(fv, istart, cf, clock) 1,19 USE DOMAIN_DECOMP, ONLY: GRID, GET, AM_I_ROOT Use GEOS_IOMod, only: GETFILE, Free_file, GEOS_VarWrite, Write_parallel USE RESOLUTION, only: IM, JM, LM, LS1 Use MODEL_COM, only: sige, sig, Ptop, DT, PMTOP Use MODEL_COM, only: U, V, T, P, PSFMPT, Q USE GEOM, only: AREAG, DXYP Use Constant, only: omega, radius, grav, rgas, kapa, deltx Type (FV_Core), Intent(InOut) :: fv integer, intent(in) :: istart Type (ESMF_Config), Intent(InOut) :: cf Type (ESMF_Clock), Intent(In) :: clock Character(Len=ESMF_MAXSTR) :: rst_file Integer :: unit Integer, Parameter :: N_TRACERS = 0 Integer :: j_0, j_1, j_0h, j_1h, L Logical :: exist Real*8 :: ak(size(sige)), bk(size(sige)) real*8, allocatable, dimension(:,:,:) :: U_d real*8, allocatable, dimension(:,:,:) :: V_d real*8, allocatable, dimension(:,:,:) :: PE, PKZ, PT ! 1) Create layout resource file - independent of actual restart file. If (AM_I_ROOT()) Then Call Write_Layout(FVCORE_LAYOUT, fv) End If Call ESMF_ConfigGetAttribute(cf, value=rst_file, label='FVCORE_INTERNAL_RESTART_FILE:', & & default=FVCORE_INTERNAL_RESTART,rc=rc) if(istart .ge. extend_run) then ! Check to see if restart file exists inquire(file=FVCORE_INTERNAL_RESTART, EXIST=exist) if (exist) then if (AM_I_ROOT()) then print*,'Using checkpoint file: ',FVCORE_INTERNAL_RESTART end if return end if ! Uh oh call stop_model('fv part of restart file not found',255) return end if ! If we got to here, then this means then we'll have to create a restart file ! from scratch. unit = GetFile(rst_file, form="unformatted", rc=rc) VERIFY_(rc) ! 1) Start date Call write_start_date(clock, unit) ! 2) Grid size Call WRITE_PARALLEL( (/ IM, JM, LM, LS1-1, N_TRACERS /), unit ) ! 3) Pressure coordinates ! Keep in mind that L is reversed between these two models Call Compute_ak_bk(ak, bk, sige, Ptop, PSFMPT, unit) Call WRITE_PARALLEL( ak, unit) Call WRITE_PARALLEL( bk, unit) Call GET(grid, j_strt=j_0, j_stop=j_1, j_strt_halo=j_0h, j_stop_halo=j_1h) ! 4) 3D fields velocities Allocate(U_d(IM, J_0:J_1, LM)) Allocate(V_d(IM, J_0:J_1, LM)) write(*,*)'Calling ComputeRestartVelocities()' Call ComputeRestartVelocities(unit, grid, U, V, U_d, V_d) !!$ call set_zonal_flow(U_d, V_d, j_0, j_1) Call GEOS_VarWrite(unit, grid % ESMF_GRID, U_d(:,J_0:J_1,:)) Call GEOS_VarWrite(unit, grid % ESMF_GRID, V_d(:,J_0:J_1,:)) Deallocate(V_d) Deallocate(U_d) ! Compute potential temperature from modelE (1 mb -> 1 pa ref) Allocate(PT(IM, J_0:J_1, LM)) Call ConvertPotTemp_GISS2FV(VirtualTemp(T(:,J_0:J_1,:), Q(:,J_0:J_1,:)), PT) Call GEOS_VarWrite(unit, grid % ESMF_GRID, PT) Deallocate(PT) ! Compute PE, PKZ from modelE Allocate(PKZ(IM, J_0:J_1, LM)) Allocate(PE(IM, J_0:J_1, LM+1)) Call ComputePressureLevels(unit, grid, VirtualTemp(T, Q), P, SIG, SIGE, Ptop, KAPA, PE, PKZ ) Call GEOS_VarWrite(unit, grid % ESMF_GRID, PE) Call GEOS_VarWrite(unit, grid % ESMF_GRID, PKZ) Deallocate(PE) Deallocate(PKZ) Call Free_File(unit) CONTAINS ! Computes virtual pot. temp. from pot. temp. and specific humidity !------------------------------------------------------------------ Function VirtualTemp(T, Q) Result(T_virt) 2,1 Use Constant, only: DELTX Real*8, Intent(In) :: T(:,:,:) Real*8, Intent(In) :: Q(:,:,:) Real*8 :: T_virt(size(T,1), size(T,2), size(T,3)) T_virt = T * (1 + deltx * Q) End Function VirtualTemp Subroutine ComputePressureLevels(unit, grid, T_virt, P, sig, sige, ptop, kapa, PE, PKZ) 1,5 USE DOMAIN_DECOMP, only: dist_grid, get USE RESOLUTION, only: IM, LM, LS1 Integer, intent(in) :: unit type (dist_grid) :: grid real*8, dimension(:,grid % j_strt_halo:,:) :: T_virt real*8, dimension(:,grid % j_strt_halo:) :: P real*8, dimension(IM,grid % j_strt:grid % j_stop,LM) :: PKZ real*8, dimension(IM,grid % j_strt:grid % j_stop,LM+1) :: PE real*8 :: sig(:), sige(:) real*8 :: ptop, kapa Integer :: L, L_fv Integer :: j_0, j_1, j_0h, j_1h Integer :: J Real*8, Allocatable :: PK(:,:,:), PELN(:,:,:), PE_trans(:,:,:) ! Request local bounds from modelE grid. Call GET(grid, j_strt=j_0, j_stop=j_1, j_strt_halo=j_0h, j_stop_halo=j_1h) PE = -99999 Allocate(pe_trans(IM,LM+1,j_0h:j_1h),pk(IM,j_0h:j_1h,LM+1), peln(IM,LM+1,j_0h:j_1h)) Call ConvertPressure_GISS2FV(EdgePressure_GISS(), PE(:,j_0:j_1,:)) ! Compute PKZ - Consistent mean for p^kappa ! use pkez() routine within the FV dycore component. Do j = j_0, j_1 PE_trans(:,:,j) = PE(:,j,:) END DO pk=0 peln=0 pkz=0 CALL pkez(1, IM, LM, J_0, J_1, 1, LM, 1, IM, PE_trans(:,:,j_0:j_1), & & PK(:,j_0:j_1,:), KAPA, LS1-1, PELN(:,:,j_0:j_1), pkz(:,j_0:j_1,:), .true.) deallocate(pe_trans, pk, peln) End Subroutine ComputePressureLevels Subroutine ComputeRestartVelocities(unit, grid, U_b, V_b, U_d, V_d) 1,4 use domain_decomp, only: DIST_GRID, NORTH, HALO_UPDATE Integer, intent(in) :: unit Type (Dist_Grid) :: grid real*8, dimension(:,grid % j_strt_halo:,:) :: U_b, V_b real*8, intent(out) :: U_d(:,:,:), V_d(:,:,:) Call HALO_UPDATE(grid, U_b, FROM=NORTH) Call HALO_UPDATE(grid, V_b, FROM=NORTH) Call Regrid_B_to_D(Reverse(U_b), Reverse(V_b), U_d, V_d) End Subroutine ComputeRestartVelocities Subroutine Compute_ak_bk(ak, bk, sige, Ptop, PSFMPT, unit) 1,1 USE RESOLUTION, only: LM, LS1 Real*8 :: sige(:) Real*8 :: Ptop, PSFMPT Integer :: unit Real*8, intent(out) :: ak(size(sige)), bk(size(sige)) Integer :: L, L_fv Do L = 1, LM+1 L_fv = LM+2-L Select Case(L) Case (:LS1-1) ak(L_fv) = Ptop*(1-sige(L)) * PRESSURE_UNIT_RATIO ! convert from hPa bk(L_fv) = sige(L) Case (LS1:) ak(L_fv) = (sige(L)*PSFMPT + Ptop) * PRESSURE_UNIT_RATIO! convert from hPa bk(L_fv) = 0 End Select End Do End Subroutine Compute_ak_bk Subroutine write_layout(fname, fv) 1,2 USE RESOLUTION, only: IM, JM, LM, LS1 Use MODEL_COM, only: DT character(len=*), intent(in) :: fname type (fv_core), intent(in) :: fv Type (esmf_axisindex), pointer :: AI(:,:) Integer :: unit Integer :: npes call esmf_vmget(fv % vm, petcount = npes, rc=rc) Allocate(AI(npes,3)) call esmf_gridgetallaxisindex(fv % grid, globalai=ai, horzrelloc=ESMF_CELL_CENTER, & & vertRelLoc=ESMF_CELL_CELL, rc=rc) unit = GetFile(fname, form="formatted", rc=rc) write(unit,*)' # empty line #' Write(unit,*)'xy_yz_decomp:',1,npes,npes,1 Write(unit,*)' im: ',IM Write(unit,*)' jm: ',JM Write(unit,*)' km: ',LM Write(unit,*)' dt: ',DT Write(unit,*)'nsplit: ',0 Write(unit,*)' ntotq: ',1 Write(unit,*)' nq: ',1 Write(unit,*)' imxy: ',IM Write(unit,'(a,100(1x,I4))')' jmxy: ',AI(:,2) % MAX-AI(:,2) % MIN+1 Write(unit,'(a,100(1x,I4))')' jmyz: ',AI(:,2) % MAX-AI(:,2) % MIN+1 Write(unit,*)' kmyz: ',LM Call Free_File(unit) Deallocate(AI) End Subroutine write_layout Subroutine write_start_date(clock, unit) 1,1 Type (ESMF_Clock), Intent(In) :: clock Integer , Intent(In) :: unit Type (ESMF_Time) :: currentTime Integer :: int_pack(6), YEAR, MONTH, DAY, HOUR, MINUTE, SECOND call ESMF_ClockGet(clock, currTime=currentTime, rc=rc) call ESMF_TimeGet(currentTime, YY=year, MM= month, & DD=day, H=hour, M=minute, & S=second, rc=rc) INT_PACK(1) = YEAR INT_PACK(2) = MONTH INT_PACK(3) = DAY INT_PACK(4) = HOUR INT_PACK(5) = MINUTE INT_PACK(6) = SECOND Call WRITE_PARALLEL(INT_PACK(1:6), unit=UNIT) End Subroutine write_start_date End Subroutine Create_Restart_File subroutine reset_tmom 1,4 USE MODEL_COM, only : im,jm,lm,t USE DOMAIN_DECOMP, ONLY: grid USE SOMTQ_COM, only : mz,tmom,qmom USE DYNAMICS, only : pmid,pedn implicit none integer :: i,j,l REAL*8 :: rdsig INTEGER :: I_0, I_1, J_1, J_0 !**** Extract useful local domain parameters from "grid" I_0 = grid%I_STRT I_1 = grid%I_STOP J_0 = grid%J_STRT J_1 = grid%J_STOP tmom = 0 ! except for vertical slopes: DO J=J_0,J_1 DO I=1,IM RDSIG=(PMID(1,I,J)-PEDN(2,I,J))/(PMID(1,I,J)-PMID(2,I,J)) TMOM(MZ,I,J,1)=(T(I,J,2)-T(I,J,1))*RDSIG DO L=2,LM-1 RDSIG=(PMID(L,I,J)-PEDN(L+1,I,J))/(PMID(L-1,I,J)-PMID(L+1,I,J)) TMOM(MZ,I,J,L)=(T(I,J,L+1)-T(I,J,L-1))*RDSIG END DO RDSIG=(PMID(LM,I,J)-PEDN(LM+1,I,J))/(PMID(LM-1,I,J)-PMID(LM,I,J)) TMOM(MZ,I,J,LM)=(T(I,J,LM)-T(I,J,LM-1))*RDSIG END DO END DO return end subroutine reset_tmom subroutine reset_qmom 1,4 USE MODEL_COM, only : im,jm,lm,q USE DOMAIN_DECOMP, ONLY: grid USE SOMTQ_COM, only : mz,tmom,qmom USE DYNAMICS, only : pmid,pedn implicit none integer :: i,j,l REAL*8 :: rdsig INTEGER :: I_0, I_1, J_1, J_0 !**** Extract useful local domain parameters from "grid" I_0 = grid%I_STRT I_1 = grid%I_STOP J_0 = grid%J_STRT J_1 = grid%J_STOP qmom = 0 ! except for vertical slopes: DO J=J_0,J_1 DO I=1,IM RDSIG=(PMID(1,I,J)-PEDN(2,I,J))/(PMID(1,I,J)-PMID(2,I,J)) QMOM(MZ,I,J,1)=(Q(I,J,2)-Q(I,J,1))*RDSIG IF(Q(I,J,1)+QMOM(MZ,I,J,1).LT.0.) QMOM(MZ,I,J,1)=-Q(I,J,1) DO L=2,LM-1 RDSIG=(PMID(L,I,J)-PEDN(L+1,I,J))/(PMID(L-1,I,J)-PMID(L+1,I,J)) QMOM(MZ,I,J,L)=(Q(I,J,L+1)-Q(I,J,L-1))*RDSIG IF(Q(I,J,L)+QMOM(MZ,I,J,L).LT.0.) QMOM(MZ,I,J,L)=-Q(I,J,L) END DO RDSIG=(PMID(LM,I,J)-PEDN(LM+1,I,J))/(PMID(LM-1,I,J)-PMID(LM,I,J)) QMOM(MZ,I,J,LM)=(Q(I,J,LM)-Q(I,J,LM-1))*RDSIG IF(Q(I,J,LM)+QMOM(MZ,I,J,LM).LT.0.) QMOM(MZ,I,J,LM)=-Q(I,J,LM) END DO END DO return end subroutine reset_qmom function PKZ_GISS() Result(PKZ) 2,10 USE RESOLUTION, only: IM, LM, LS1 Use MODEL_COM, only : SIG, Ptop, PSFMPT, P use DOMAIN_DECOMP, only: grid, get USE CONSTANT, only: KAPA REAL*8 :: PKZ(IM,grid % J_STRT:grid % J_STOP,LM) INTEGER :: L, j_0, j_1 call get(grid, J_STRT=J_0, J_STOP=J_1) Do L = 1, LM If (L < LS1) THEN PKZ(:,:,L) = (SIG(L)*P(:,J_0:J_1) + Ptop) ** KAPA Else PKZ(:,:,L) = (SIG(L)*PSFMPT + Ptop) ** KAPA End IF End Do end function PKZ_GISS ! Convert Potential Temperature into (dry) Temperature function DryTemp_GISS() Result(T_dry) 4,10 USE RESOLUTION, only: IM, LM, LS1 Use MODEL_COM, only: T USE DOMAIN_DECOMP, only: grid, GET REAL*8 :: T_dry(IM,grid % J_STRT:grid % J_STOP,LM) REAL*8 :: PKZ(IM,grid % J_STRT:grid % J_STOP,LM) INTEGER :: J_0,J_1 Call Get(grid, J_STRT=J_0, J_STOP=J_1) PKZ = PKZ_GISS() T_dry = PKZ * T(:,J_0:J_1,:) end function DryTemp_GISS Subroutine Copy_modelE_to_FV_import(fv) 1,5 USE MODEL_COM, only: Q ! Secific Humidity USE MODEL_COM, only: ZATMO ! Geopotential Height? USE MODEL_COM, Only : U, V, T Use DOMAIN_DECOMP, Only: grid, Get Type (FV_CORE) :: fv Integer :: nq Integer :: j_0, j_1 Call Get(grid, j_strt=j_0, j_stop=j_1) ! 1) Link/copy modelE data to import state fv % phis=ZATMO(:,j_0:j_1) #ifdef NO_FORCING fv % phis = 0 #endif ! Moisture fv % Q = Reverse(Q(:,j_0:j_1,:)) #ifdef NO_FORCING fv % Q = 0 #endif End Subroutine Copy_modelE_to_FV_import Subroutine Copy_FV_export_to_modelE(fv) 1,14 Use Resolution, only: IM,JM,LM,LS1 USE DYNAMICS, ONLY: PUA,PVA,SDA, PU, PV, SD Use MODEL_COM, only: U, V, T, P, PSFMPT, Q Use MODEL_COM, only : Ptop, P USE DOMAIN_DECOMP, only: grid, GET, SOUTH, NORTH, HALO_UPDATE Use ATMDYN, only: CALC_AMPK USE GEOM Type (FV_CORE) :: fv real*4, Dimension(:,:,:), Pointer :: U_a, V_a, T_fv, PLE Integer :: unit Integer :: i,j,k Integer :: i_0, i_1, j_0, j_1 Call Get(grid, i_strt=i_0, i_stop=i_1, j_strt=j_0, j_stop=j_1) ! First compute updated values for modelE. Then capture the ! new state in fv % *_old for computing tendencies. ! ---------------------------------------------------------- ! Velocity field !--------------- call ESMFL_StateGetPointerToData ( fv % export,U_a,'U',rc=rc) VERIFY_(rc) call ESMFL_StateGetPointerToData ( fv % export,V_a,'V',rc=rc) VERIFY_(rc) call Regrid_A_to_B(Reverse(U_a), Reverse(V_a), U(:,J_0:J_1,:), V(:,J_0:J_1,:)) fv % U_old = U(:,J_0:J_1,:) fv % V_old = V(:,J_0:J_1,:) ! Potential temperature (save dry Temperature for computing tendencies) !---------------------------------------------------------------------- call ESMFL_StateGetPointerToData ( fv % export,T_fv,'TH',rc=rc) VERIFY_(rc) call ConvertPotTemp_FV2GISS(T_fv, T) ! Use edge pressure export to update surface pressure !---------------------------------------------------------------------- call ESMFL_StateGetPointerToData ( fv % export,PLE,'PLE',rc=rc) VERIFY_(rc) call ConvertPressure_FV2GISS(PLE, fv % PE_old) ! Just need surface pressure - Ptop P(:,J_0:J_1) = fv % PE_old(:,:,1) - Ptop CALL CALC_AMPK(LS1-1) ! Preserve state information for later computation of tendencies. fv % dT_old = DryTemp_GISS() fv % dPT_old = DeltPressure_DryTemp_GISS() #if defined(USE_FV_Q) Q(:,j_0:j_1,:) = Reverse(fv % Q) #endif End Subroutine Copy_FV_export_to_modelE subroutine ConvertPressure_GISS2FV_r4(P_giss, P_fv) 1 Real*8, intent(in) :: P_giss(:,:,:) Real*4, intent(out) :: P_fv(:,:,:) ! no halo in this case P_fv = Reverse(P_giss * PRESSURE_UNIT_RATIO) end Subroutine ConvertPressure_GISS2FV_r4 subroutine ConvertPressure_GISS2FV_r8(P_giss, P_fv) 1 Real*8, intent(in) :: P_giss(:,:,:) Real*8, intent(out) :: P_fv(:,:,:) ! no halo in this case P_fv = Reverse(P_giss * PRESSURE_UNIT_RATIO) end Subroutine ConvertPressure_GISS2FV_r8 ! Convert pressure from GISS representation to FV representation. ! Both the order of levels and the units must be adjusted. subroutine ConvertPressure_FV2GISS(P_fv, P_giss) 1 Real*4, intent(in) :: P_fv(:,:,:) Real*8, intent(out) :: P_giss(:,:,:) P_giss = Reverse(P_fv / PRESSURE_UNIT_RATIO) end Subroutine ConvertPressure_FV2GISS ! Convert potential temperature between the two representations. subroutine CnvPotTemp_GISS2FV_r8(PT_giss, PT_fv) 1,1 Use CONSTANT, only : KAPA Real*8, intent(in) :: PT_giss(:,:,:) Real*8, intent(out) :: PT_fv(:,:,:) PT_fv = Reverse(PT_giss * REF_RATIO ** KAPA) end Subroutine CnvPotTemp_GISS2FV_r8 ! Convert potential temperature between the two representations. subroutine CnvPotTemp_GISS2FV_r4(PT_giss, PT_fv) 1,1 Use CONSTANT, only : KAPA Real*8, intent(in) :: PT_giss(:,:,:) Real*4, intent(out) :: PT_fv(:,:,:) PT_fv = Reverse(PT_giss * REF_RATIO ** KAPA) end Subroutine CnvPotTemp_GISS2FV_r4 !------------------------------------------------------------------------ ! Convert potential temperature as exported by FV to the corresponding ! potential temperauture within modelE. Note that FV export uses ! a reference pressure of 10^5 pa. !------------------------------------------------------------------------ subroutine CnvPotTemp_FV2GISS_r4(PT_fv, PT_giss) 1,1 Use CONSTANT, only : KAPA Real*4, intent(in) :: PT_fv(:,:,:) Real*8, intent(out) :: PT_giss(:,:,:) ! As an export, FV provides Pot Temp with a reference ! pressure of 10^5 Pa, whereas GISS uses 100 Pa (1 mb) REAL*8, PARAMETER :: REF_RATIO = 100000 / 100 integer :: n ! PT_GISS has a halo, while PT_fv does not. n = size(PT_GISS,2) ! exclude halo PT_giss(:,2:n-1,:) = Reverse(PT_fv / REF_RATIO ** KAPA) end Subroutine CnvPotTemp_FV2GISS_r4 !------------------------------------------------------------------------ ! Given velocities from modelE (B grid, k=1 bottom), produce ! velocities for import to FV (A grid, k=1 top) !------------------------------------------------------------------------ subroutine ConvertUV_GISS2FV(U_giss, V_giss, U_fv, V_fv) 1,1 Real*8, intent(in) :: U_giss(:,:,:) Real*8, intent(in) :: V_giss(:,:,:) Real*4, intent(out) :: U_fv(:,:,:) Real*4, intent(out) :: V_fv(:,:,:) Call Regrid_B_to_A(Reverse(U_giss), Reverse(V_giss), U_fv, V_fv) end subroutine ConvertUV_GISS2FV !------------------------------------------------------------------------ ! This following routine interpolates U, V from the Arakawa A grid ! to the Arakawa B grid. B-grid velocities correspond to GISS modelE, ! while A-grid velocities are used for import/export of the FV component. !------------------------------------------------------------------------ subroutine Regrid_A_to_B(U_a, V_a, U_b, V_b) 1,5 Use Resolution, only : IM, LM Use DOMAIN_DECOMP, only : grid, get, SOUTH, HALO_UPDATE Real*4, intent(in), Dimension(:,grid % J_STRT:,:) :: U_a, V_a Real*8, intent(out),Dimension(:,grid % J_STRT:,:) :: U_b, V_b Real*8, allocatable, dimension(:,:,:) :: Ua_halo, Va_halo Integer :: i,j,k,im1 integer :: j_0stgr, j_1stgr,J_0h,J_1h,J_0,J_1 Call Get(grid, J_STRT_STGR=j_0stgr, J_STOP_STGR=j_1stgr, J_STRT_HALO=J_0H, J_STOP_HALO=J_1H, & & J_STRT=J_0, J_STOP=J_1) Allocate(Ua_halo(IM,J_0h:J_1h,LM), Va_halo(IM,J_0h:J_1h,LM)) Ua_halo(:,J_0:J_1,:) = U_a Va_halo(:,J_0:J_1,:) = V_a Call Halo_Update(grid, Ua_halo,FROM=SOUTH) Call Halo_Update(grid, Va_halo,FROM=SOUTH) Do k = 1, LM Do j = j_0STGR, j_1STGR im1 = IM Do i = 1, IM u_b(i,j,k) = (Ua_halo(im1,j-1,k) + Ua_halo(i,j-1,k) + Ua_halo(im1,j,k) + Ua_halo(i,j,k))/4 v_b(i,j,k) = (Va_halo(im1,j-1,k) + Va_halo(i,j-1,k) + Va_halo(im1,j,k) + Va_halo(i,j,k))/4 im1 = i End do end do end do Deallocate(Ua_halo, Va_halo) end subroutine Regrid_A_to_B !------------------------------------------------------------------------ ! This following routine interpolates U, V from the Arakawa B grid ! to the Arakawa A grid. B-grid velocities correspond to GISS modelE, ! while A-grid velocities are used for import/export of the FV component. !------------------------------------------------------------------------ subroutine Regrid_B_to_A(U_b, V_b, U_a, V_a) 1,7 Use Resolution, only : IM, JM, LM Use DOMAIN_DECOMP, only : grid, get, NORTH, HALO_UPDATE Real*8, intent(in), Dimension(:,grid % J_STRT:,:) :: U_b, V_b Real*4, intent(out), Dimension(:,grid % J_STRT:,:) :: U_a, V_a Real*8, allocatable, dimension(:,:,:) :: Ub_halo, Vb_halo Integer :: i,j,k,im1 integer :: j_0, j_1, j_0s, j_1s, j_0h, j_1h logical :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE Call Get(grid, J_STRT=J_0, J_STOP=J_1, J_STRT_SKP=J_0s, J_STOP_SKP=J_1S, & & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, HAVE_NORTH_POLE=HAVE_NORTH_POLE, & & J_STRT_HALO=J_0H, J_STOP_HALO=J_1H) Allocate(Ub_halo(IM,J_0h:J_1h,LM), Vb_halo(IM,J_0h:J_1h,LM)) Ub_halo(:,J_0:J_1,:) = U_b Vb_halo(:,J_0:J_1,:) = V_b Call Halo_Update(grid, Ub_halo,FROM=NORTH) Call Halo_Update(grid, Vb_halo,FROM=NORTH) Do k = 1, LM Do j = j_0s,j_1s im1 = IM Do i = 1, IM u_a(im1,j,k) = (ub_halo(im1,j,k) + ub_halo(i,j,k) + ub_halo(im1,j+1,k) + ub_halo(i,j+1,k))/4 v_a(im1,j,k) = (vb_halo(im1,j,k) + vb_halo(i,j,k) + vb_halo(im1,j+1,k) + vb_halo(i,j+1,k))/4 im1 = i End do end do end do deallocate(ub_halo, Vb_halo) ! Polar conditions are a bit more complicated ! First determine an "absolute" U, V at the pole, then use sin/cos to ! map to the indidual longitudes. If (HAVE_SOUTH_POLE) then Call FixPole(U_b(:,2,:), V_b(:,2,:), U_a(:,1,:), V_a(:,1,:)) End If If (HAVE_NORTH_POLE) then Call FixPole(U_b(:,JM,:), V_b(:,JM,:), U_a(:,JM,:), V_a(:,JM,:)) End If Contains ! This routine works for both poles. ! The sign is correct for the NP, but (-1) appears as a factor twice ! in the derivation for the south pole. !---------------------------------------- Subroutine FixPole(Ub, Vb, Ua, Va) 4,2 USE GEOM, only : SINIV, COSIV, SINIP, COSIP ! trig of lon and stgr lon Real*8, intent(in) :: Ub(IM, LM), Vb(IM,LM) Real*4, intent(out) :: Ua(IM, LM), Va(IM,LM) Real*8 :: us, vs Integer :: k do k = 1, LM us = Sum( -ub(:,k) * SINIP - vb(:,k) * COSIP) / IM vs = Sum( ub(:,k) * COSIP - vb(:,k) * SINIP) / IM ua(:,k) = -us * SINIV + vs * COSIV va(:,k) = -us * COSIV - vs * SINIV end do End Subroutine FixPole end subroutine Regrid_B_to_A !-------------------------------------------------------------------- ! This following routine interpolates U, V from the Arakawa B grid ! to the Arakawa D grid. B-grid velocities correspond to GISS modelE, ! while D-grid velocities are used for the initial conditions of FV !-------------------------------------------------------------------- Subroutine regrid_B_to_D(u_b, v_b, u_d, v_d) 1,5 USE RESOLUTION, only: IM, JM, LM Use DOMAIN_DECOMP, only: grid, get Implicit None Real*8, intent(in) :: U_b(:,grid % j_strt_halo:,:) Real*8, intent(in) :: V_b(:,grid % j_strt_halo:,:) Real*8, intent(out) :: U_d(:,grid % j_strt:,:) Real*8, intent(out) :: V_d(:,grid % j_strt:,:) Integer :: i, j, k, ip1 Integer :: j_0, j_1, j_0s, j_1s Logical :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE Call Get(grid, j_strt=j_0, j_stop=j_1, j_strt_skp=j_0s, j_stop_skp=j_1s, & & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, HAVE_NORTH_POLE=HAVE_NORTH_POLE) If (HAVE_SOUTH_POLE) Then ! 1st interpolate to 'A' grid Call FixPole(U_b(:,2,:), V_b(:,2,:), U_d(:,1,:)) ! V_d(:,1,:) = (V_b(:,2,:) + CSHIFT(V_b(:,2,:),1,1))/2 U_d(:,1,:)=0 ! not used (but needs legal value) End If If (HAVE_NORTH_POLE) Then ! 1st interpolate to 'A' grid Call FixPole(U_b(:,JM,:), V_b(:,JM,:), U_d(:,JM,:)) ! V_d(:,JM,:) = (V_b(:,JM,:) + CSHIFT(V_b(:,JM,:),1,1))/2 End If ! V-grid Do k = 1, LM Do j = j_0s, j_1s Do i = 1, IM V_d(i,j,k) = (V_b(i,j,k) + V_b(i,j+1,k)) /2 End Do ! V_d(:,j,:) = (V_b(:,j,:) + CSHIFT(V_b(:,j,:),1,1))/2 End Do End Do ! U-grid Do k = 1, LM Do j = j_0s, j_1 i = IM Do ip1 = 1, IM U_d(i,j,k) = (U_b(i,j,k) + U_b(ip1,j,k)) /2 i = ip1 End Do ! U_d(:,j,:) = (U_b(:,j,:) + CSHIFT(U_b(:,j,:),1,1))/2 End Do End Do Contains Subroutine FixPole(Ub, Vb, Ud) 4,2 USE GEOM, only : SINIP, COSIP ! trig of lon and stgr lon Real*8, intent(in) :: Ub(IM, LM), Vb(IM,LM) Real*8, intent(out) :: Ud(IM, LM) Real*8 :: US, VS Integer :: k do k = 1, LM us = Sum( -ub(:,k) * SINIP - vb(:,k) * COSIP) / IM vs = Sum( ub(:,k) * COSIP - vb(:,k) * SINIP) / IM ud(:,k) = (ub(:,k) + (-us * SINIP + vs * COSIP))/2 end do End Subroutine FixPole End Subroutine regrid_B_to_D ! Reverse the order of vertical levels. Function reverse_3d_r8(A) Result(B) 1 Real*8, Intent(In) :: A(:,:,:) Real*8 :: B(Size(A,1),Size(A,2),size(A,3)) Integer, parameter :: K_IDX = 3 Integer :: k, n n = Size(A, K_IDX) Do k = 1, n B(:,:,k) = A(:,:,1+n-k) End Do End Function reverse_3d_r8 ! Single precision variant of Reverse() Function reverse_3d_r4(A) Result(B) 1 Real*4, Intent(In) :: A(:,:,:) Real*4 :: B(Size(A,1),Size(A,2),size(A,3)) Integer, parameter :: K_IDX = 3 Integer :: k, n n = Size(A, K_IDX) Do k = 1, n B(:,:,k) = A(:,:,1+n-k) End Do End Function reverse_3d_r4 Subroutine Write_Profile(arr, name),6 Use RESOLUTION, Only: IM, JM, LM Use DOMAIN_DECOMP, Only: grid, PACK_DATA, AM_I_ROOT Real*8, intent(in) :: arr(:,:,:) character(len=*), intent(in) :: name Integer :: k, km Real*8 :: rng(3,LM) Real*8 :: arr_global(IM,JM,size(arr,3)) Real*8 :: arr_tmp(IM,grid % j_STRT_HALO:grid % J_stop_HALO,size(arr,3)) arr_tmp(:,grid % J_strt:grid % J_STOP,:)=arr Call PACK_DATA(grid, arr_tmp, arr_global) IF (AM_I_ROOT()) Then rng(1,:) = MINVAL(MINVAL(arr_global,DIM=1),DIM=1) rng(2,:) = MAXVAL(MAXVAL(arr_global,DIM=1),DIM=1) rng(3,:) = SUM(SUM(arr_global,DIM=1),DIM=1)/(IM*JM) print*,'***********' print*,'stats for ',trim(name) km = size(arr,3) Do k = 1, km Write(*,'(a,i4.0,3(f21.9,1x))')'k:',k,rng(:,k) End Do print*,'***********' print*,' ' End IF End Subroutine Write_Profile function compute_phi(P, T, SZ, zatmo) result(phi) 1,6 USE CONSTANT, only: KAPA, BYKAPA, BYKAPAP1, BYKAPAP2, RGAS USE MODEL_COM, only: LS1, LM, DSIG, SIG, SIGE, PTOP, PSFMPT USE MODEL_COM, only: IM, JM, LM USE DOMAIN_DECOMP, Only: grid, Get USE GEOM, only: IMAXJ implicit none real*8, intent(in) :: P(:,grid % J_STRT_HALO:) real*8, intent(in) :: T(:,grid % J_STRT_HALO:,:) real*8, intent(in) :: SZ(:,grid % J_STRT_HALO:,:) real*8, intent(in) :: ZATMO(:,grid % J_STRT_HALO:) real*8 :: phi(IM,grid % J_STRT_HALO:grid % J_STOP_HALO,LM) REAL*8 :: PKE(LS1:LM+1) REAL*8 :: PIJ, PHIDN REAL*8 :: DP, BYDP, P0, TZBYDP, X REAL*8 :: PUP, PKUP, PKPUP, PKPPUP REAL*8 :: PDN, PKDN, PKPDN, PKPPDN INTEGER :: I, J, L INTEGER :: J_0, J_1 LOGICAL :: HAVE_NORTH_POLE, HAVE_SOUTH_POLE call get(grid, J_STRT=J_0, J_STOP=J_1, HAVE_NORTH_POLE=HAVE_NORTH_POLE, & & HAVE_SOUTH_POLE = HAVE_SOUTH_POLE) DO L=LS1,LM+1 PKE(L)=(PSFMPT*SIGE(L)+PTOP)**KAPA END DO !$OMP PARALLEL DO PRIVATE(I,J,L,DP,P0,PIJ,PHIDN,TZBYDP,X, !$OMP* BYDP,PDN,PKDN,PKPDN,PKPPDN,PUP,PKUP,PKPUP,PKPPUP) DO J=J_0,J_1 DO I=1,IMAXJ(J) PIJ=P(I,J) PDN=PIJ+PTOP PKDN=PDN**KAPA PHIDN=ZATMO(I,J) !**** LOOP OVER THE LAYERS DO L=1,LM PKPDN=PKDN*PDN PKPPDN=PKPDN*PDN IF(L.GE.LS1) THEN DP=DSIG(L)*PSFMPT BYDP=1./DP P0=SIG(L)*PSFMPT+PTOP TZBYDP=2.*SZ(I,J,L)*BYDP X=T(I,J,L)+TZBYDP*P0 PUP=SIGE(L+1)*PSFMPT+PTOP PKUP=PKE(L+1) PKPUP=PKUP*PUP PKPPUP=PKPUP*PUP ELSE DP=DSIG(L)*PIJ BYDP=1./DP P0=SIG(L)*PIJ+PTOP TZBYDP=2.*SZ(I,J,L)*BYDP X=T(I,J,L)+TZBYDP*P0 PUP=SIGE(L+1)*PIJ+PTOP PKUP=PUP**KAPA PKPUP=PKUP*PUP PKPPUP=PKPUP*PUP END IF !**** CALCULATE PHI, MASS WEIGHTED THROUGHOUT THE LAYER PHI(I,J,L)=PHIDN+RGAS*(X*PKDN*BYKAPA-TZBYDP*PKPDN*BYKAPAP1 & & -(X*(PKPDN-PKPUP)*BYKAPA-TZBYDP*(PKPPDN-PKPPUP)*BYKAPAP2) & & *BYDP*BYKAPAP1) !**** CALULATE PHI AT LAYER TOP (EQUAL TO BOTTOM OF NEXT LAYER) PHIDN=PHIDN+RGAS*(X*(PKDN-PKUP)*BYKAPA-TZBYDP*(PKPDN-PKPUP) & & *BYKAPAP1) PDN=PUP PKDN=PKUP END DO END DO END DO !$OMP END PARALLEL DO !**** SET POLAR VALUES FROM THOSE AT I=1 IF (HAVE_SOUTH_POLE) THEN DO L=1,LM PHI(2:IM,1,L)=PHI(1,1,L) END DO END IF IF (HAVE_NORTH_POLE) THEN DO L=1,LM PHI(2:IM,JM,L)=PHI(1,JM,L) END DO END IF end function compute_phi subroutine clear_accumulated_mass_fluxes() 1,1 USE DYNAMICS, ONLY: PUA,PVA,SDA PUA(:,:,:) = 0.0 PVA(:,:,:) = 0.0 SDA(:,:,:) = 0.0 end subroutine clear_accumulated_mass_fluxes subroutine accumulate_mass_fluxes(fv) 1,7 Use Resolution, only: IM,JM,LM,LS1 USE GEOM, ONLY: DYP, DXV, DXYP, IMAXJ, BYIM USE DYNAMICS, ONLY: PUA,PVA,SDA USE DYNAMICS, ONLY: PU,PV,CONV,SD,PIT USE MODEL_COM, only: DTsrc,DT,DSIG USE DOMAIN_DECOMP, only: get, grid, NORTH, SOUTH, HALO_UPDATE ! Use Constant, only: radius,pi,grav implicit none type (FV_core) :: fv real*4, Dimension(:,:,:), Pointer :: PLE, mfx_X, mfx_Y, mfx_Z integer :: J_0, J_1 integer :: J_0S, J_1S integer :: J_0H, J_1H integer :: i,im1,j,l,k integer :: rc logical :: HAVE_NORTH_POLE, HAVE_SOUTH_POLE real*8 :: DTLF, DTfac real*8 :: area,dlon,dlat,acap,rcap real*8 :: rX,rY,rZ real*8 :: sum1, sum2, mySum real*8, allocatable :: sine(:),cosp(:),cose(:) REAL*8 PVS,PVN REAL*8 PVSA(LM),PVNA(LM) REAL*8 TMPim1(IM) real*8, parameter :: grav=9.80 real*8, parameter :: pi=3.14159265358979323846 real*8, parameter :: radius= 6376000.00000000 DTfac = DT Call Get(grid, j_strt=j_0, j_stop=j_1, J_STRT_SKP=J_0S, J_STOP_SKP=J_1S, & & J_STRT_HALO=J_0H, J_STOP_HALO = J_1H, & & HAVE_NORTH_POLE = HAVE_NORTH_POLE, HAVE_SOUTH_POLE = HAVE_SOUTH_POLE) ! Horizontal and Vertical mass fluxes !--------------- ! \item {\tt MFX}: Mass-Weighted U-Wind on C-Grid (Pa m^2/s) ! \item {\tt MFY}: Mass-Weighted V-wind on C-Grid (Pa m^2/s) ! \item {\tt MFZ}: Vertical mass flux (kg/(m^2*s)) call ESMFL_StateGetPointerToData ( fv % export,mfx_X,'MFX',rc=rc) VERIFY_(rc) call ESMFL_StateGetPointerToData ( fv % export,mfx_Y,'MFY',rc=rc) VERIFY_(rc) call ESMFL_StateGetPointerToData ( fv % export,mfx_Z,'MFZ',rc=rc) VERIFY_(rc) call ESMFL_StateGetPointerToData ( fv % export,PLE,'PLE',rc=rc) VERIFY_(rc) allocate ( sine(jm) ) allocate ( cosp(jm) ) allocate ( cose(jm) ) dlon = 2.0*pi/im dlat = pi/(jm-1) do j=2,jm sine(j) = sin(-0.5*pi + ((j-1)-0.5)*(pi/(jm-1))) enddo cosp( 1) = 0.0 cosp(jm) = 0.0 do j=2,jm-1 cosp(j) = (sine(j+1)-sine(j)) / dlat enddo do j=2,jm cose(j) = 0.5 * (cosp(j-1) + cosp(j)) enddo cose(1) = cose(2) PU(1:IM,J_0:J_1,1:LM) = mfx_X PV(1:IM,J_0:J_1,1:LM) = mfx_Y #ifdef NO_MASS_FLUX mfx_X = 0 mfx_Y = 0 mfx_Z = 0 #endif PU = Reverse(PU)/PRESSURE_UNIT_RATIO PV = Reverse(PV)/PRESSURE_UNIT_RATIO mfx_Z = Reverse(mfx_Z)/PRESSURE_UNIT_RATIO ! Adjust area scale factors between FV and GISS do l=1,lm do j=j_0,j_1 do i=1,im PU(i,j,l) = PU(i,j,l)*DYP(j)/(radius*dlat) enddo do i=1,im PV(i,j,l) = PV(i,j,l)*DXV(j)/(cose(j)*radius*dlon) enddo enddo enddo ! Shift C-grid PU to Eastward orientation for GISS do l=1,lm do j=j_0,j_1 im1 = 1 do i=im,1,-1 TMPim1(i) = PU(im1,j,l) im1 = i enddo PU(:,j,l) = TMPim1 enddo enddo ! Change Units of vertical mass fluxes do l=0,lm do j=j_0,j_1 area = DXYP(j) do i=1,im mfx_Z(i,j-j_0+1,l) = grav*area*mfx_Z(i,j-j_0+1,l) ! convert to (mb m^2/s) enddo enddo enddo SD(1:IM,J_0:J_1,1:LM-1) = (mfx_Z(:,:,1:LM-1)) ! SD only goes up to LM-1 ! Surface Pressure tendency - vert integral of horizontal convergence PIT(1:IM,J_0:J_1) = mfx_Z(:,:,0) + sum(SD(1:IM,J_0:J_1,1:LM-1),3) deallocate ( sine ) deallocate ( cosp ) deallocate ( cose ) ! Recopy into CONV to support prior usage CONV(:,J_0:J_1,1) = PIT(:,J_0:J_1) CONV(:,J_0:J_1,2:LM) = SD(:,J_0:J_1,1:LM-1) PUA(:,J_0:J_1,:) = PUA(:,J_0:J_1,:) + PU(:,J_0:J_1,:)*DTfac PVA(:,J_0:J_1,:) = PVA(:,J_0:J_1,:) + PV(:,J_0:J_1,:)*DTfac SDA(:,J_0:J_1,1:LM-1) = SDA(:,J_0:J_1,1:LM-1) + SD(:,J_0:J_1,1:LM-1)*DTfac end subroutine accumulate_mass_fluxes subroutine set_zonal_flow(U_d, V_d, j0, j1),2 use GEOM, only: cosp use DOMAIN_DECOMP, only: grid Real*8, intent(out) :: U_d(:,grid % j_strt:,:) Real*8, intent(out) :: V_d(:,grid % j_strt:,:) integer :: j0, j1 integer :: j do j = j0, j1 U_d(:,j,:) = cosp(j) V_d(:,j,:) = 0 end do end subroutine set_zonal_flow subroutine printSnapshot(fv),1 use MODEL_COM, only: U, V, T, Q, P type (fv_core) :: fv integer, save :: counter = 0 counter = counter + 1 write(40,*)'Iteration: ',counter write(40,*)'U : ', U(4, 3, 3) write(40,*)'V : ', V(4, 3, 3) write(40,*)'T : ', T(4, 3, 3) write(40,*)'Q : ', Q(4, 3, 3) write(40,*)'P : ', P(4, 3) write(40,*)'dudt : ', fv % dudt(4, 3, 3), sum(fv% dudt) write(40,*)'dvdt : ', fv % dvdt(4, 3, 3), sum(fv% dvdt) write(40,*)'dtdt : ', fv % dtdt(4, 3, 3), sum(fv% dtdt) write(40,*)'dpedt : ', fv % dpedt(4, 3, 14), sum(fv% dpedt) write(40,*)'Q : ', fv % q(4, 3, 3) write(40,*)'PHIS : ', fv % phis(4, 3) write(40,*)'*******' write(40,*)' ' end subroutine printSnapshot end module FV_INTERFACE_MOD !---------------------------------------------------------------- ! The following routine interfaces the VERIFY_ macro with the ! GISS termination routine. The line number and return code are ! written into a buffer which is passed to stop_model(). !---------------------------------------------------------------- Subroutine abort_core(line,rc) 1,1 Implicit None Integer, Intent(In) :: line Integer, Intent(In) :: rc Character(len=100) :: buf Write(buf,*)'FV core failure at line',line,'\n error code:', rc Call stop_model(Trim(buf),99) End Subroutine abort_core module dynamics_save implicit none private public :: initialize public :: save_dynamics_state public :: restore_dynamics_state public :: write_dynamics_state real*8, allocatable, dimension(:,:,:), save :: U_save, V_save, T_save real*8, allocatable, dimension(:,:), save :: P_save real*8, allocatable, dimension(:,:,:), save :: Q_save, WM_save real*8, allocatable, dimension(:,:,:), save :: gz_save, phi_save real*8, allocatable, dimension(:,:,:,:), save :: tmom_save real*8, allocatable, dimension(:,:,:), save :: pdsig_save contains subroutine initialize() 1,2 use resolution, only : IM, JM, LM use somtq_com, only : tmom, NMOM allocate(U_save(IM,1:JM,LM)) allocate(V_save(IM,1:JM,LM)) allocate(T_save(IM,1:JM,LM)) allocate(P_save(IM,1:JM)) allocate(Q_save(IM,1:JM,LM)) allocate(WM_save(IM,1:JM,LM)) allocate(gz_save(IM,1:JM,LM)) allocate(phi_save(IM,1:JM,LM)) allocate(tmom_save(NMOM,IM,1:JM,LM)) allocate(pdsig_save(LM,IM,1:JM)) end subroutine initialize subroutine save_dynamics_state(),4 use resolution, only : IM, JM, LM use model_com, only: U, V, T, P, Q, WM use dynamics, only: gz, phi, pdsig use somtq_com, only : tmom U_save = U(:,1:JM,:) V_save = V(:,1:JM,:) T_save = T(:,1:JM,:) P_save = P(:,1:JM) Q_save = Q(:,1:JM,:) WM_save = WM(:,1:JM,:) gz_save = gz(:,1:JM,:) phi_save = phi(:,1:JM,:) tmom_save = tmom(:,:,1:JM,:) pdsig_save = pdsig(:,:,1:JM) end subroutine save_dynamics_state subroutine restore_dynamics_state(),4 use resolution, only : IM, JM, LM use model_com, only: U, V, T, P, Q, WM use dynamics, only: gz, phi, pdsig use somtq_com, only : tmom U(:,1:JM,:) = U_save V(:,1:JM,:) = V_save T(:,1:JM,:) = T_save P(:,1:JM) = P_save Q(:,1:JM,:) = Q_save WM(:,1:JM,:) = WM_save gz(:,1:JM,:) = gz_save phi(:,1:JM,:) = phi_save tmom(:,:,1:JM,:) = tmom_save pdsig(:,:,1:JM) = pdsig_save end subroutine restore_dynamics_state subroutine write_dynamics_state(dyn),9 use resolution, only : IM, JM, LM use model_com, only: U, V, T, P character(len=*), intent(in) :: dyn integer :: unit select case(dyn) case ('FV') unit = 21 case ('MODELE A') unit = 22 case ('MODELE B') unit = 23 case ('MODELE C') unit = 24 case default end select write(unit) 'U',U(:,1:JM,:) write(unit) 'V',V(:,1:JM,:) write(unit) 'T',T(:,1:JM,:) write(unit) 'P',P(:,1:JM) select case(dyn) case ('FV') unit = 31 case ('MODELE A') unit = 32 case ('MODELE B') unit = 33 case ('MODELE C') unit = 34 case default end select write(*,*) 'DYNAMICS for ' // trim(dyn) call write_surface_zone('U',U(:,:,1)) call write_surface_zone('V',V(:,:,1)) call write_surface_zone('T',T(:,:,1)) call write_surface_zone('P',P(:,:)) call write_avg('U',U) call write_avg('V',V) call write_avg('T',T) contains subroutine write_surface_zone(name, array) 4 character(len=*), intent(in) :: name real*8, intent(in) :: array(:,:) integer :: j write(unit,*)' ' write(unit,*)'*********************************' write(unit,*)'Surface Zonal Average for ' // trim(name) write(unit,*)'*********************************' do j = 1, JM write(unit,*) j, sum(array(:,j))/IM end do write(unit,*)' ' end subroutine write_surface_zone subroutine write_avg(name, array) 3,1 use geom, only: AREAG, dxyp character(len=*), intent(in) :: name real*8, intent(in) :: array(:,:,:) integer :: k write(unit,*)' ' write(unit,*)'*********************************' write(unit,*)' Level Average for ' // trim(name) write(unit,*)'*********************************' do k = 1, LM write(unit,*) k, sum(sum(array(:,1:JM,k),1)*dxyp(1:jm))/AREAG end do write(unit,*)' ' end subroutine write_avg end subroutine write_dynamics_state end module dynamics_save #else module FV_INTERFACE_MOD 8,10 implicit none private public :: DryTemp_GISS, Write_Profile contains Subroutine Write_Profile(arr, name),6 Use RESOLUTION, Only: IM, JM, LM Use DOMAIN_DECOMP, Only: grid, PACK_DATA, AM_I_ROOT Real*8, intent(in) :: arr(:,:,:) character(len=*), intent(in) :: name Integer :: k, km Real*8 :: rng(3,LM) Real*8 :: arr_global(IM,JM,size(arr,3)) Real*8 :: arr_tmp(IM,grid % j_STRT_HALO:grid % J_stop_HALO,size(arr,3)) arr_tmp(:,grid % J_strt:grid % J_STOP,:)=arr Call PACK_DATA(grid, arr_tmp, arr_global) IF (AM_I_ROOT()) Then rng(1,:) = MINVAL(MINVAL(arr_global,DIM=1),DIM=1) rng(2,:) = MAXVAL(MAXVAL(arr_global,DIM=1),DIM=1) rng(3,:) = SUM(SUM(arr_global,DIM=1),DIM=1)/(IM*JM) print*,'***********' print*,'stats for ',trim(name) km = size(arr,3) Do k = 1, km Write(*,'(a,i4.0,3(f21.9,1x))')'k:',k,rng(:,k) End Do print*,'***********' print*,' ' End IF End Subroutine Write_Profile function PKZ_GISS() Result(PKZ) 2,10 USE RESOLUTION, only: IM, LM, LS1 Use MODEL_COM, only : SIG, Ptop, PSFMPT, P use DOMAIN_DECOMP, only: grid, get USE CONSTANT, only: KAPA REAL*8 :: PKZ(IM,grid % J_STRT:grid % J_STOP,LM) INTEGER :: L, j_0, j_1 call get(grid, J_STRT=J_0, J_STOP=J_1) Do L = 1, LM If (L < LS1) THEN PKZ(:,:,L) = (SIG(L)*P(:,J_0:J_1) + Ptop) ** KAPA Else PKZ(:,:,L) = (SIG(L)*PSFMPT + Ptop) ** KAPA End IF End Do end function PKZ_GISS ! Convert Potential Temperature into (dry) Temperature function DryTemp_GISS() Result(T_dry) 4,10 USE RESOLUTION, only: IM, LM, LS1 Use MODEL_COM, only: T USE DOMAIN_DECOMP, only: grid, GET REAL*8 :: T_dry(IM,grid % J_STRT:grid % J_STOP,LM) REAL*8 :: PKZ(IM,grid % J_STRT:grid % J_STOP,LM) INTEGER :: J_0,J_1 Call Get(grid, J_STRT=J_0, J_STOP=J_1) PKZ = PKZ_GISS() T_dry = PKZ * T(:,J_0:J_1,:) end function DryTemp_GISS end module FV_INTERFACE_MOD #endif