! $Id: POS_GEOS5PlugMod.F90,v 1.3.2.1 2007/05/24 18:29:12 pschopf Exp $ #include "MAPL_Generic.h" ! GEOS-5 default real kind module POS_GEOS5PlugMod ! A MAPL/ESMF Gridded Component tha acts as a wrapper for Poseidon. ! It uses ESMF AND MAPL. It has heavy dependencies on GMAO_gems and Poseidon. use ESMF_Mod use MAPL_Mod ! Used in Initialize: use Communication_Primitives, only: GEMS_cp_initarch => cp_Initarch use GHOST_MODULE, only:& HaloFill=>GHOST,& ! GHOST Communication Routine gmNORTH=>NORTH,& ! Constant indicating "NORTH" gmEAST=>EAST,& ! Constant indicating "EAST" gmWEST=>WEST,& ! Constant indicating "WEST" gmSOUTH=>SOUTH,& ! Constant indicating "SOUTH" gmA_STAGGER=>mass_point,& gmB_STAGGER=>z_point,& gmU_STAGGER=>u_point,& gmV_STAGGER=>v_point use global_constants, only: CTOKELV use Poseidon_Grid_Module, only: & T_Poseidon_Grid, & POSEIDON_Init_grid_gridfile => Init_grid_gridfile, & POSEIDON_Init_grid_restart => Init_grid_restart, & PGRID_GetWorldLatLons, & Reset_grid_hmins, & Destroy_grid, & PGRID_findcore, & pg_C_GRID, & pg_GRIDCENTER, & pg_SURFACE use Poseidon_Restart_Module, only: & Initialize_poseidon => Init,& Print_Run_Parameters, & Save_restart use Poseidon_Module, only: & T_poseidon_ocean => T_ocean,& T_poseidon_forcing, & Run_poseidon => run, & Create_poseidon_forcing => Create_forcing, & MarkDiagForHistory, & EnableDiagnostics, & DisableDiagnostics, & ResetDiagnostic, & Destroy_ocean, & Destroy_forcing, & Validate_ocean, & PrintWanted use Poseidon_History_Module, only: & T_poseidon_history,& Create_ohist,& Destroy_ohist,& Initialize_ohist,& Write_ohist,& Accumulate_ohist use Zeus_clocks, only: & Z_clock => T_clock,& Z_alarm => T_alarm,& Z_DateTime => T_DateTime, & Z_TimeInterval => T_TimeInterval, & Z_DateTimeInterval => T_DateTimeInterval, & Z_as_DateTime => as_DateTime, & Z_as_TimeInterval => as_TimeInterval, & Z_as_Seconds => as_Seconds, & Z_print_clock, & operator(+), & operator(-), & Z_Init_Alarm, & Z_Set_Interval,& Z_Set_Time,& Z_Set_Alarm_Ringing, & Z_tick, & Z_Free_Clock,& Z_Free_Alarm,& Z_Shut_Alarm,& Z_Ringing use EQ_STATE_MODULE, only: & dbdtemp, & dbdsalt, & RHO_OCN_R, & RHOCP_R, & CP_OCN #ifdef ASSIM_MODE use odas_iau_module, only: t_odas_iau, create_odas_iau, apply_odas_iau #endif ! Nothing on the Poseidon side is visible through this module. implicit none private ! Only public things are the IRF routines. public :: SetServices ! These are the Poseidon objects type POS_MAPL_Type type (T_poseidon_ocean ) :: ocean type (T_poseidon_grid ) :: grid type (T_poseidon_forcing) :: oforce type (T_poseidon_history) :: ohist type (Z_clock) :: zclock Type (Z_Alarm ) :: EOR_Alarm ! Alarm for stopping Poseidon Type (Z_Alarm ) :: history_Alarm ! Alarm for writing history Type (Z_Alarm ) :: history_Acc_Alarm ! Alarm for accumulating history averages integer :: history_unit ! Fortran unit number for writing history data integer :: grads_descriptor_unit ! Fortran unit number for writing GrADS control logical :: write_control_file = .true. ! Should we write GrADS control? type (Z_DateTimeInterval) :: historyInterval character(ESMF_MAXSTR) :: ohist_template #ifdef ASSIM_MODE type (t_odas_iau) :: iau ! Incremental analysis updating object #endif end type POS_MAPL_Type type POS_MAPLWrap_Type type(POS_MAPL_Type), pointer :: Ptr end type POS_MAPLWrap_Type ! ErrLog Variables character(len=ESMF_MAXSTR) :: IAm integer :: STATUS character(len=ESMF_MAXSTR) :: COMP_NAME interface UpdateFromConfig module procedure UpdateRealFromConfig module procedure UpdateIntFromConfig module procedure UpdateLogFromConfig end interface UpdateFromConfig interface Assignment(=) module procedure Z_to_E_TimeInterval module procedure E_to_Z_TimeInterval module procedure E_to_Z_DateTime module procedure Z_to_E_DateTime end interface integer, parameter :: realHIGH = SELECTED_REAL_KIND(12,307) ! kind of old REAL*8 contains !BOP ! !IROUTINE: SetServices -- Sets ESMF services for for POS_Plug wrapper ! !INTERFACE: subroutine SetServices ( GC, RC ) ! !ARGUMENTS: type(ESMF_GridComp), intent(INOUT) :: GC ! gridded component integer, optional :: RC ! return code ! !DESCRIPTION: The SetServices for the Poseidon\_Plug needs to register its ! Initialize, Run and Finalize. It has no children. !EOP !============================================================================= ! ! Locals !============================================================================= ! 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 !BOP ! !IMPORT STATE: call MAPL_AddImportSpec(GC, & SHORT_NAME = 'TAUX', & LONG_NAME = 'Agrid_eastward_stress_on_skin', & UNITS = 'N m-2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'TAUY', & LONG_NAME = 'Agrid_northward_stress_on_skin', & UNITS = 'N m-2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'USTAR3', & LONG_NAME = 'friction velocity cubed', & UNITS = 'm+3 s-3', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'PS', & LONG_NAME = 'Surface Atmospheric Pressure', & UNITS = 'Pa', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'SWHEAT', & LONG_NAME = 'solar_heating_rate', & UNITS = 'W m-2', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'QFLX', & LONG_NAME = 'freshwater_flux_from_skin_to_ocean',& UNITS = 'kg m-2 s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'HFLX', & LONG_NAME = 'total_heat_flux_into_ocean', & UNITS = 'W m-2', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'SFLX', & LONG_NAME = 'salt_flux_into_ocean', & UNITS = 'kg m-2 s-1', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'TR', & LONG_NAME = 'tracer_mixing_ratios', & UNITS = '1', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & STAT = MAPL_BundleItem, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddImportSpec(GC, & SHORT_NAME = 'TRFLUX', & LONG_NAME = 'surface_fluxes_of_tracers_into_ocean',& UNITS = 'X', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & STAT = MAPL_BundleItem, & RC=STATUS ) VERIFY_(STATUS) ! !EXPORT STATE: call MAPL_AddExportSpec(GC, & SHORT_NAME = 'US', & LONG_NAME = 'top_layer_Agrid_eastward_velocity', & UNITS = 'm s-1 ', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'VS', & LONG_NAME = 'top_layer_Agrid_northward_velocity',& UNITS = 'm s-1 ', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'TS', & LONG_NAME = 'top_layer_temperature', & UNITS = 'K', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'SS', & LONG_NAME = 'top_layer_salinity', & UNITS = 'psu', & DIMS = MAPL_DimsHorzOnly, & VLOCATION = MAPL_VLocationNone, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'DH', & LONG_NAME = 'layer_mass', & UNITS = 'dyn m', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) call MAPL_AddExportSpec(GC, & SHORT_NAME = 'MASK', & LONG_NAME = '3D_land-sea mask', & UNITS = 'none', & DIMS = MAPL_DimsHorzVert, & VLOCATION = MAPL_VLocationCenter, & RC=STATUS ) VERIFY_(STATUS) !EOP ! Set the Initialize, Run, Finalize entry points ! ---------------------------------------------- call MAPL_GridCompSetEntryPoint ( GC, ESMF_SETINIT, Initialize, status) VERIFY_(STATUS) call MAPL_GridCompSetEntryPoint ( GC, ESMF_SETRUN, Run, status) VERIFY_(STATUS) call MAPL_GridCompSetEntryPoint ( GC, ESMF_SETFINAL, Finalize, status) VERIFY_(STATUS) ! Add the Profiling timers ! ------------------------ call MAPL_TimerAdd(GC, name="INITIALIZE" ,RC=STATUS) VERIFY_(STATUS) call MAPL_TimerAdd(GC, name="RUN" ,RC=STATUS) VERIFY_(STATUS) call MAPL_TimerAdd(GC, name="FINALIZE" ,RC=STATUS) VERIFY_(STATUS) call MAPL_TimerAdd(GC, name="--PosRun" ,RC=STATUS) VERIFY_(STATUS) call MAPL_TimerAdd(GC, name="--PosHistory",RC=STATUS) VERIFY_(STATUS) ! Generic SetServices -- call Last ! ------------------- call MAPL_GenericSetServices ( GC, RC=STATUS ) VERIFY_(STATUS) ! All done ! -------- RETURN_(ESMF_SUCCESS) end subroutine SetServices !BOP ! !IROUTINE: Initialize -- Initialize method for the for POS_Plug wrapper ! !INTERFACE: subroutine Initialize ( GC, IMPORT, EXPORT, CLOCK, RC ) ! !ARGUMENTS: type(ESMF_GridComp), intent(INOUT) :: GC ! Gridded component type(ESMF_State), intent(INOUT) :: IMPORT ! Import state type(ESMF_State), intent(INOUT) :: EXPORT ! Export state type(ESMF_Clock), intent(INOUT) :: CLOCK ! The clock integer, optional, intent( OUT) :: RC ! Error code: ! !DESCRIPTION: ! The Initialize method for Poseidon\_Plug performs the ! following tasks. ! \begin{enumerate} ! \item MAPL\_Generic initialization. This must come first. ! \item Instantiates the POS\_MAPL\_internal\_state ! \item Initialize GMAO\_gems communications. This needs to ! match the layout, virtual machine, and MPI Communicator from ! the ESMF\_Grid assigned to this gridded component. ! \item Initializes the Poseidon grid in the GMAO\_gems communicator. ! The Poseidon grid information can come from a variety of sources - ! the restart or a special grid file. The ESMF Configuration ! file is used to pass this information. ! \item Initializes the Poseidon state from its private restart data. At ! the present, only the ZDF data format is permitted for restart ! initialization. The name of the restart file is obtrined from the ! Configuration file. ! \item Synchronizes the Poseidon clock with the ESMF clock, sets alarms, etc. ! This initialization assumes that the current time on the ESMF clock ! represents the initial instant from which the run will be made. ! \item Sets up diagnostics and private Poseidon-style history. Poseidon ! diagnostics are extensive and go well beyond what is available ! in the state. At the present time, they can only be computed ! and saved using the Poseidon private history mechanism. A mechanism ! to export any poseidon diagnostic so that MAPL History can report it ! is under development. ! \item Validates the Poseidon state and parameters as suitable for running. ! \item Sets export fields based on the restart conditions. Things like ! sea surface temperature need to be up-to-date. ! \end{enumerate} ! ! Because Poseidon was designed as an object-based model, it has ! always required a supervisory driver program to perform similar ! initialization functions. By making the POS\_Plug wrapper, the ! Initialize routine is now required to perform these same initializations. ! There are numerous internal routines which were lifted from earlier ! supervisory codes to perform the real work in Initialize. !EOP ! Locals with ESMF and MAPL types type(ESMF_DELayout) :: Layout type(ESMF_VM) :: VM type(MAPL_MetaComp), pointer :: MAPL type(ESMF_Config) :: CF type(ESMF_Grid) :: ESMFGrid type(ESMF_Time) :: MyTime ! Locals type (T_poseidon_ocean ), pointer :: ocean type (T_poseidon_grid ), pointer :: grid type (T_poseidon_forcing), pointer :: oforce type (T_poseidon_history), pointer :: ohist type (Z_clock), pointer :: zclock ! Exports - We need to fill these at startup real, pointer :: TS (:,:) real, pointer :: SS (:,:) real, pointer :: DH (:,:,:) real, pointer :: MASK(:,:,:) type(POS_MAPL_Type), pointer :: POS_MAPL_internal_state type(POS_MAPLWrap_Type) :: wrap character(ESMF_MAXSTR) :: OceanGridFile ! Name of the poseidon grid file character(ESMF_MAXSTR) :: OceanRestartFile ! Name of the poseidon restart file character(ESMF_MAXSTR) :: OceanRestartType ! Type of the poseidon restart (zdf) integer :: DT_OCEAN integer :: dimcount integer :: comm logical :: has_extmode integer :: YEAR,MONTH,DAY,HR,MN,SC integer :: i1,i2,j1,j2 ! Begin... ! Get the target components name and set-up traceback handle. ! ----------------------------------------------------------- Iam = "Poseidon_Initialize" call ESMF_GridCompGet( gc, NAME=comp_name, RC=status ) VERIFY_(STATUS) Iam = trim(comp_name) // trim(Iam) ! Get my internal MAPL_Generic state (Only used for profilers in this routine) !----------------------------------- call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS) VERIFY_(STATUS) ! Profilers !---------- call MAPL_TimerOn(MAPL,"TOTAL" ) call MAPL_TimerOn(MAPL,"INITIALIZE") ! Generic initialize ! ------------------ call MAPL_GenericInitialize( GC, IMPORT, EXPORT, CLOCK, RC=status ) VERIFY_(STATUS) ! Get the grid, configuration !---------------------------- call ESMF_GridCompGet( GC, grid=ESMFGrid, CONFIG = CF, RC=status ) VERIFY_(STATUS) ! Get the layout from the grid !----------------------------- call ESMF_GridGet(ESMFGrid, deLayout=Layout,rc=STATUS) VERIFY_(STATUS) ! Get the dimensions of the DElayout !----------------------------------- call ESMF_DELayoutGet(Layout, dimCount=dimcount, RC=STATUS ) VERIFY_(STATUS) ASSERT_(dimcount == 2) call ESMF_DELayoutGetVm(Layout, VM=VM, RC=STATUS ) VERIFY_(STATUS) ! Set the time for Poseidon !--------------------- call ESMF_ClockGet(CLOCK, currTIME=MyTime, RC=STATUS) VERIFY_(STATUS) call ESMF_TimeGet (MyTime, & YY=YEAR, MM=MONTH, DD=DAY, & H=HR, M =MN, S =SC, & RC=STATUS ) VERIFY_(STATUS) ! Allocate this instance of the internal state and put it in wrapper. ! ------------------------------------------------------------------- allocate( POS_MAPL_internal_state, stat=status ) VERIFY_(STATUS) wrap%ptr => POS_MAPL_internal_state ! Save pointer to the wrapped internal state in the GC ! ---------------------------------------------------- call ESMF_UserCompSetInternalState ( GC, 'POS_MAPL_state',wrap,status ) VERIFY_(STATUS) ocean => POS_MAPL_internal_state%ocean oforce => POS_MAPL_internal_state%oforce ohist => POS_MAPL_internal_state%ohist grid => POS_MAPL_internal_state%grid zclock => POS_MAPL_internal_state%zclock ! GEMS initialization !------------------------------------------------------ call ESMF_VMGet(VM, mpiCommunicator=Comm, rc=STATUS) VERIFY_(STATUS) call GEMS_cp_Initarch(Comm) ! Do Poseidon stuff !------------- call ESMF_ConfigGetAttribute( CF, OceanGridFile, & LABEL=trim(comp_name)//"_GRID_FILE:", & RC=status) VERIFY_(status) call write_parallel( 'OceanGridFile: ['//trim(OceanGridFile)//']' ) call ESMF_ConfigGetAttribute( CF, OceanRestartType, & LABEL=trim(comp_name)//"_POSEIDON_RESTART_TYPE:", & DEFAULT='zdf',RC=status) VERIFY_(status) call ESMF_ConfigGetAttribute( CF, OceanRestartFile,& LABEL=trim(comp_name)//"_POSEIDON_RESTART_FILE:", & DEFAULT='notFound',RC=status) VERIFY_(status) ! Initialize DE layout, a poseidon grid, and an ESMF grid: ! (everything needed to get ourselves going on a distributed grid) call Initialize_grids_and_layouts( RC=status ) VERIFY_(status) i1 = grid%i1 i2 = grid%i2 j1 = grid%j1 j2 = grid%j2 call Initialize_ocean( rc=status ) VERIFY_(status) call Setup_clocks( rc=status ) VERIFY_(status) call Choose_diagnostics( rc=status ) VERIFY_(status) call Setup_history( rc=status ) VERIFY_(status) call Z_PRINT_CLOCK( zclock ) call Validate_ocean( ocean, QUIET=.false.) ! Profilers ! --------- call MAPL_TimerOff(MAPL,"INITIALIZE") call MAPL_TimerOff(MAPL,"TOTAL" ) call MAPL_GetPointer(EXPORT, TS, 'TS' , alloc=.true., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, SS, 'SS' , alloc=.true., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, DH, 'DH' , alloc=.true., RC=STATUS) VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, MASK, 'MASK', alloc=.true., RC=STATUS) VERIFY_(STATUS) if(associated(TS )) then TS = ocean%state%T(i1:i2,j1:j2,1) + CTOKELV end if if(associated(SS )) then SS = ocean%state%S(i1:i2,j1:j2,1) + ocean%state%sbias end if if(associated(DH)) then DH = ocean%state%H(i1:i2,j1:j2,:) end if if(associated(MASK)) then MASK = grid%bh3(i1:i2,j1:j2,:) end if ! All Done !--------- RETURN_(ESMF_SUCCESS) contains subroutine Initialize_grids_and_layouts(RC) integer, intent(out), optional :: RC ! return code type(esmf_field) :: my_esmf_field; integer :: DIMS(2) integer, pointer :: i_per_DE_column(:) ! number of i grid-points in each DE column integer, pointer :: j_per_DE_row(:) ! number of j grid-points in each DE row integer :: nDE_cols, nDE_rows integer :: j integer :: imw,jmw,kmw ! world grid dimensions integer, allocatable :: cellCountPerDEPerDim(:,:) integer, allocatable :: cellCountPerDEcolumn(:) integer, allocatable :: cellCountPerDErow(:) call ESMF_DELayoutGet(layout, deCountPerDim = DIMS, rc=status) VERIFY_(status) nDE_cols = DIMS(1) nDE_rows = DIMS(2) allocate( i_per_DE_column(nDE_cols), j_per_DE_row(nDE_rows), STAT=status ) VERIFY_(status) allocate( cellCountPerDEPerDim( nDE_cols * nDE_rows, 2 ), STAT=status ) VERIFY_(status) call ESMF_GridGet( ESMFGrid , ESMF_CELL_CENTER, cellCountPerDEPerDim=cellCountPerDEPerDim, rc=status) VERIFY_(status) allocate( cellCountPerDEcolumn(nDE_cols), cellCountPerDErow(nDE_rows), STAT=status) VERIFY_(status) cellCountPerDEcolumn(1:nDE_cols) = cellCountPerDEPerDim(1:nDE_cols,1) do j=1,nDE_rows cellCountPerDErow(j) = cellCountPerDEPerDim( j*nDE_cols-1,2 ) enddo call Write_parallel(trim(comp_name) // ' ESMF Grid Cell Count Layout:') call Write_parallel(cellCountPerDEcolumn(1:nDE_cols),format='(" IMS: ",(256I4))') call Write_parallel(cellCountPerDErow(1:nDE_rows),format='(" JMS: ",(256I4))') call Write_parallel(' ') deallocate( cellCountPerDEPerDim, stat=status) VERIFY_(status) ! Look at the gridfile and extract im,jm,km, then ! consult the config file and get the DE layout information. ! Determine the size of the grid call Get_Dimensions(imw,jmw,kmw,i_per_DE_column,j_per_DE_row,RC=status) VERIFY_(status) ASSERT_( imw == sum(cellCountPerDEcolumn) ) ASSERT_( jmw == sum(cellCountPerDErow) ) if ( .not. all(cellCountPerDEcolumn(1:nDE_cols) == i_per_DE_column(1:nDE_cols)) ) then call write_parallel('NOTE: i_per_DE_column being changed to match ESMF grid') i_per_DE_column(:) = cellCountPerDEcolumn(:) endif if ( .not. all( cellCountPerDErow(:) == j_per_DE_row(:)) ) then call write_parallel('NOTE: j_per_DE_row being changed to match ESMF grid') j_per_DE_row(:) = cellCountPerDErow(:) endif !----------------------------------- if ( trim(OceanGridFile) == 'none' ) then if ( trim(OceanRestartFile) .ne. 'notFound' ) then if ( OceanRestartType == 'zdf' .or. OceanRestartType == 'ZDF' ) then call POSEIDON_Init_grid_restart(& trim(OceanRestartFile),& grid , & EXTMODE=has_extmode,& IMS=i_per_DE_column, & JMS=j_per_DE_row ) else call write_parallel(' Request to get Poseidon grid from restart but restart is not .zdf') ASSERT_(.false.) endif else ! Here we would look to get a grid from a new restart file ASSERT_(.false.) endif else call POSEIDON_Init_grid_gridfile( OceanGridFile, kmw, & grid , & EXTMODE=has_extmode,& DEFAULT_LAYOUT = pg_C_GRID, & IMS=i_per_DE_column, & JMS=j_per_DE_row ) endif ! Now check that the ESMF grid agrees with the poseidon grid. call Check_Grids( grid, ESMFgrid, RC=status) VERIFY_(status) deallocate( i_per_DE_column, j_per_DE_row, stat=status) VERIFY_(status) deallocate( cellCountPerDEcolumn, cellCountPerDErow, stat=status) VERIFY_(status) RETURN_(status) end subroutine Initialize_grids_and_layouts ! Read the ocean grid size from the header !----------------------------------------- subroutine Get_Dimensions( imw, jmw, kmw, i_per_DE_column,j_per_DE_row, RC ) integer, intent(out) :: imw ! Size of world grid inner dimension integer, intent(out) :: jmw ! Size of world grid outer dimension integer, intent(out) :: kmw ! Number of layers integer, pointer :: i_per_DE_column(:) ! number of i grid-points in each DE column integer, pointer :: j_per_DE_row(:) ! number of j grid-points in each DE row integer, optional, intent(out) :: RC ! Read the ocean grid file header to obtain im,jm,km ! and to tell whether we are to run an external mode integer :: unit integer :: i real :: head(20) real :: x0,y0,dx,dy integer :: nDE_cols, nDE_rows if ( trim(OceanGridFile) .ne. 'none' ) then call write_parallel( 'Getting dimensions from gridfile: ['//trim(OceanGridFile)//']' ) unit=Getfile( trim(OceanGridFile), form='unformatted', RC=status) VERIFY_(status) call Read_parallel( layout, head, unit=unit, RC=status) VERIFY_(status) call Free_file(unit) ! once the header is read, go ahead and close the unit imw = head(9)-1 ! The first record on the gridfile has an extra row jmw = head(10)-1 ! and column (it holds corner information) kmw = head(11) x0 = head(12) y0 = head(13) dx = head(14) dy = head(15) else call write_parallel( 'Getting dimensions from restartfile: ['//trim(OceanRestartFile)//']' ) unit=Getfile( trim(OceanRestartFile), form='unformatted', RC=status) VERIFY_(status) call Read_parallel( layout, head, unit=unit, RC=status) VERIFY_(status) call Free_file(unit) ! once the header is read, go ahead and close the unit imw = head(9) ! The first record on the gridfile has an extra row jmw = head(10) ! and column (it holds corner information) kmw = head(11) dx = head(14) dy = head(15) x0 = head(12) - dx ! poseidon.restarts.f90 puts the NE corner into these header vars y0 = head(13) - dy ! we need to re-cast to the SW corner endif ! Set up the decomposition information. We either get the data ! from the config file, or use the GEOS method for doing this nDE_cols = ubound( i_per_DE_column, 1 ) call ESMF_ConfigFindLabel( cf, trim(comp_name)//'_I_PER_DE_COL:', RC = status ) if( STATUS == ESMF_SUCCESS ) then VERIFY_(status) do i = 1, nDE_cols call ESMF_ConfigGetAttribute( cf, i_per_DE_column(i), RC = STATUS ) VERIFY_(status) enddo else i_per_DE_column(:) = 0 call set_ij_per_proc( imw, i_per_DE_column, nDE_cols ) endif ASSERT_( imw == sum(i_per_DE_column) ) call Write_parallel(trim(comp_name) // ' PE Layout:') call Write_parallel(i_per_DE_column(1:nDE_cols),format='(" IMS: ",(256I4))') call Write_parallel(' ') ! Get the JMS vector ! ------------------ nDE_rows = ubound( j_per_DE_row, 1 ) call ESMF_ConfigFindLabel( cf, trim(comp_name)//'_J_PER_DE_ROW:', RC = status ) if( status == ESMF_SUCCESS ) then do i = 1, nDE_rows call ESMF_ConfigGetAttribute( cf, j_per_DE_row(i), RC = status ) VERIFY_(status) enddo else j_per_DE_row(:) = 0 call set_ij_per_proc( jmw, j_per_DE_row, nDE_rows ) endif ASSERT_( jmw == sum(j_per_DE_row) ) call Write_parallel(trim(comp_name) // ' PE Layout:') call Write_parallel(j_per_DE_row(1:nDE_rows),format='(" JMS: ",(256I4))') call Write_parallel(' ') call ESMF_ConfigFindLabel( cf, trim(comp_name)//'_EXTMODE:', RC = status ) if( status == ESMF_SUCCESS ) then call ESMF_ConfigGetAttribute( cf, i, RC = status ) has_extmode = i > 0 else has_extmode = .true. endif if ( present(RC) ) RC=ESMF_SUCCESS end subroutine Get_Dimensions ! Check that the world grid points assigned to each processor ! is consistent, or if they are all zero, then come up with ! a reasonable distribution subroutine set_ij_per_proc( npoints, points_per_proc, nprocs ) integer, intent(in) :: npoints integer, intent(inout) :: points_per_proc(:) integer, intent(in) :: nprocs integer :: boundaries(0:SIZE(points_per_proc)),i if ( SUM(points_per_proc(1:nprocs) ) > 0 ) then ASSERT_( SUM(points_per_proc(1:nprocs)) == npoints ) else boundaries = (/(i,i=0,nprocs)/) * REAL(npoints / nprocs) boundaries(nprocs) = npoints points_per_proc(1:nprocs) = boundaries(1:nprocs)-boundaries(0:nprocs-1) endif end subroutine set_ij_per_proc subroutine Check_Grids( pgrid, ESMFgrid, rc) use GLOBAL_CONSTANTS, only: DTOR type(T_Poseidon_grid), intent(IN ) :: pgrid ! Poseidon grid object type (ESMF_Grid), intent(IN ) :: ESMFgrid integer, intent( OUT) :: rc ! Local vars integer :: status character(len=ESMF_MAXSTR), parameter :: IAm='Poseidon_Initialize.Check_grids' integer :: L character(len=ESMF_MAXSTR) :: gridname real, allocatable :: lats(:), lons(:) real(ESMF_KIND_R8) :: deltaZ real, pointer :: xne_egrid(:,:), yne_egrid(:,:) real, allocatable :: test(:,:),xy(:,:) real :: epsilon integer :: i,j,point(2) integer :: i1,i2,j1,j2,ie,je,ip,jp,ni,nj type (ESMF_logical) :: is_Periodic(2) real(ESMF_KIND_R8) :: twopi real :: pi pi = acos(-1.0) i1=pgrid%i1 i2=pgrid%i2 j1=pgrid%j1 j2=pgrid%j2 ni=i2-i1+1 nj=j2-j1+1 ! Let's be sure that the grid has been initialized ASSERT_( pgrid%initialized ) call GridCornerGet(ESMFgrid, xne_egrid, & NAME="Longitude", & rc=status) VERIFY_(status) ASSERT_( lbound(xne_egrid,1) == 0 ) ASSERT_( ubound(xne_egrid,1) == ni ) ASSERT_( lbound(xne_egrid,2) == 0 ) ASSERT_( ubound(xne_egrid,2) == nj ) call WRITE_PARALLEL('POSEIDON: Verified that GEMS and ESMF grids have same size in x') allocate( test(0:ni,0:nj), xy(0:ni,0:nj), stat=status) VERIFY_(status) ! Normalize longitudes by Pi (pgrid is in degrees, ESMF_grid ! can be anything) twopi = GridTwoPi(ESMFgrid, rc=status ) VERIFY_(status) xy(0:ni,0:nj) = pgrid%lon_ne(i1-1:i2,j1-1:j2)/180. xne_egrid = xne_egrid * 2. / twopi ! to handle negative values and 2Pi wraps, add 3Pi to the ! difference, MOD with 2Pi, then subtract Pi. ! then take absolute value test = xne_egrid - xy + 3. test = abs( mod( test, 2. ) - 1. ) ! Be sure that the grid agrees to within some tolerance ! use 1% of the nominal x grid spacing epsilon = ( 0.01 / pgrid%imw ) if ( maxval(test) > epsilon ) then print *,'Error in Poseidon Initialize: ESMF grid does not look like Poseidon grid' point = maxloc(test) ie=point(1) - 1 je=point(2) - 1 print *,' i= ',ie,' j= ',je,' x(ESMF_grid):',xne_egrid(ie,je), & ' x(Pgrid): ',xy(ie,je),& ' diff: ',test(ie,je),& ' epsilon: ',epsilon,& ' twopi: ',twopi ASSERT_(.false.) else call WRITE_PARALLEL('POSEIDON: Verified that GEMS and ESMF grids have x corners') endif call GridCornerGet(ESMFgrid, yne_egrid, & NAME="Latitude", & rc=status) VERIFY_(status) ASSERT_( lbound(yne_egrid,1) == 0 ) ASSERT_( ubound(yne_egrid,1) == ni ) ASSERT_( lbound(yne_egrid,2) == 0 ) ASSERT_( ubound(yne_egrid,2) == nj ) call WRITE_PARALLEL('POSEIDON: Verified that GEMS and ESMF grids have same size in y') ! normalize latitude by Pi xy(0:ni,0:nj) = pgrid%lat_ne(i1-1:i2,j1-1:j2) / 180. yne_egrid = yne_egrid * 2. / twopi test = abs( yne_egrid - xy ) ! Be sure that the grid agrees to within some tolerance ! use 1% of the nominal x grid spacing epsilon = ( 180. / pgrid%jmw ) * 0.01 * DTOR if ( maxval(test) > epsilon ) then print *,'Error in Poseidon Initialize: ESMF grid does not look like Poseidon grid' point = maxloc(test) ie = point(1) - 1 je = point(2) - 1 print *,' i= ',ie,' j= ',je,' y(ESMF_grid):',yne_egrid(ie,je), & ' y(Pgrid): ',xy(ie,je),& ' diff: ',test(ie,je),& ' epsilon: ',epsilon ASSERT_(.false.) else call WRITE_PARALLEL('POSEIDON: Verified that GEMS and ESMF grids have y corners') endif deallocate( xne_egrid, yne_egrid, test, xy, stat=status ) VERIFY_(status) ! We should do some checking here, to be sure that the ! ESMF grid agrees with the poseidon grid, no? RETURN_(status) end subroutine Check_Grids function GridTwoPi( grid, rc ) result(twopi) type(ESMF_Grid), intent(IN ) :: GRID integer, optional :: rc real(ESMF_KIND_R8) :: twopi type(ESMF_CoordOrder) :: coordorder character(len=ESMF_MAXSTR):: IAm="Poseidon_Initialize.ESMFL_GridTwoPi" integer :: index real(ESMF_KIND_R8), dimension(ESMF_MAXGRIDDIM) :: minGlobalCoordPerDim real(ESMF_KIND_R8), dimension(ESMF_MAXGRIDDIM) :: maxGlobalCoordPerDim twopi = 0.0 call ESMF_GridGet( grid, coordorder=coordorder, & minGlobalCoordPerDim=minGlobalCoordPerDim,& maxGlobalCoordPerDim=maxGlobalCoordPerDim,& rc=status ) VERIFY_(status) if ( coordorder == ESMF_COORD_ORDER_XYZ .or. & coordorder == ESMF_COORD_ORDER_XYZ ) then index = 1 elseif ( coordorder == ESMF_COORD_ORDER_YXZ .or. & coordorder == ESMF_COORD_ORDER_ZXY ) then index = 2 elseif ( coordorder == ESMF_COORD_ORDER_YZX .or. & coordorder == ESMF_COORD_ORDER_ZYX ) then index = 3 else ASSERT_(.false.) endif twopi = maxGlobalCoordPerDim(index)-minGlobalCoordPerDim(index) RETURN_(status) end function GridTwoPi subroutine GridCornerGet(GRID, coord, name, rc) type(ESMF_Grid), intent(IN ) :: GRID real, dimension(:,:), pointer :: coord character (len=*) , intent(IN) :: name integer, optional :: rc ! local variables type (ESMF_Array), pointer :: earray(:) integer :: rank type(ESMF_DataType) :: type type(ESMF_DataKind) :: kind integer :: counts(ESMF_MAXDIM) integer :: crdOrder real(ESMF_KIND_R8), pointer :: r8d2(:,:,:) integer :: STATUS integer :: j character(len=ESMF_MAXSTR):: IAm="Poseidon_Initialize.GridCornerGet" type (ESMF_Array), target :: tarray(2) earray=>tarray call ESMF_GridGetCoord(GRID, horzRelloc=ESMF_CELL_CENTER, & cornerCoord=earray, rc=status) VERIFY_(STATUS) if (name == "Longitude") then crdOrder = 1 else if (name == "Latitude") then crdOrder = 2 else STATUS=ESMF_FAILURE VERIFY_(STATUS) endif call ESMF_ArrayGet(earray(crdOrder), RANK=rank, TYPE=type, KIND=kind, & COUNTS=counts, RC=status) VERIFY_(STATUS) if (rank /=3 ) then ! ALT: for now STATUS=ESMF_FAILURE VERIFY_(STATUS) endif if (type /= ESMF_DATA_REAL) then STATUS=ESMF_FAILURE VERIFY_(STATUS) endif if ( counts(1) /= 4 ) then !PSS: the first dimension is the 4 corners STATUS=ESMF_FAILURE VERIFY_(STATUS) endif ! allocate coord with (0:im,0:jm), so they look like the NE corners allocate( coord(0:counts(2), 0:counts(3)), STAT=status) VERIFY_(status) if (kind == ESMF_R8) then ! We are counting on ESMF returning corners numbered ! clockwise from the SW, so 1=SW, 2=SE, 3=NE, 4=NW ! There should be a way to check this. call ESMF_ArrayGetData(earray(crdOrder), r8d2, RC=status) VERIFY_(STATUS) coord(1:,1:) = r8d2(3,:,:) coord(0,0) = r8d2(1,1,1) coord(1:,0) = r8d2(2,:,1) coord(0,1:) = r8d2(4,1,:) else !!! ALT: This here is an ugly kludge. Replace ASAP coord = 0.0 ! if (crdOrder == 1) ! STATUS=ESMF_FAILURE ! VERIFY_(STATUS) ! End of kludge endif RETURN_(ESMF_SUCCESS) end subroutine GridCornerGet subroutine Initialize_ocean(rc) integer, optional, intent(out) :: rc character(ESMF_MAXSTR) :: OceanNamelistFile ! Now initialize the poseidon state ! call write_parallel('Initializing - type = '//trim(OceanRestartType) ) if ( OceanRestartType == 'zdf' .or. OceanRestartType == 'ZDF' ) then call ESMF_ConfigGetAttribute( cf, OceanNamelistFile, & LABEL=trim(comp_name)//"_POSEIDON_NAMELIST_FILE:", & DEFAULT='none',RC=status) VERIFY_(status) call Initialize_poseidon( ocean, grid, zclock, & control_file=trim(OceanNamelistFile),& restart_file=trim(OceanRestartFile),& grid_file=trim(OceanGridFile),& ID=trim(comp_name),& set_clock=.true.,& verbose=.false.) else ! In this case, the internal state was saved by GEOS ASSERT_( .false. ) !call relate_poseidon_to_internalstate( & ! POSEIDONSTATE=ocean, & ! GEOSINTERNAL=internalState, & ! RC=status) !VERIFY_(status) endif call Create_poseidon_forcing( oforce , grid, ocean%param%ntrace ) ! call WRITE_PARALLEL('Getting modified run parameters from config file') call Update_Parameters_from_Config( cf, ocean, rc=status ) VERIFY_(status) call Write_parallel('Summary of Run Parameters for Poseidon') call Print_Run_Parameters( ocean ) call Z_Set_Interval( zclock, ocean%param%delt, status=status) VERIFY_(status) end subroutine Initialize_ocean subroutine Update_Parameters_from_Config(cf, ocean, rc) type (ESMF_Config) :: cf type (T_poseidon_ocean) :: ocean integer, optional, intent(out) :: rc character(len=ESMF_MAXSTR), parameter :: IAm='Initialize_ocean.Update_Parameters_from_Config' integer :: status real :: binit(ocean%g%kmg) integer :: n call UpdateFromConfig( cf, ocean%param%delt, "TIMESTEP:", rc=status) VERIFY_(status) call Z_Set_Interval(ocean%alarm%step, ocean%param%delt, status=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%extmode, "EXTERNAL_MODE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Ext%substeps, "EXTERNAL_SUBSTEPS:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Mass_damping_time, "MASS_DAMPING_TIME:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%beta_h, "BETA_H:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%U_timescale, "U_FILTER_TIMESCALE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%U_frequency, "U_FILTER_FREQUENCY:", rc=status) VERIFY_(status) call Z_Set_Interval(ocean%alarm%ufilter, ocean%param%filt%U_frequency, status=status) ! Dont verify since status=-1 is fine here VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%U_order, "U_FILTER_ORDER:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%H_timescale, "H_FILTER_TIMESCALE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%H_frequency, "H_FILTER_FREQUENCY:", rc=status) VERIFY_(status) call Z_Set_Interval(ocean%alarm%hfilter, ocean%param%filt%H_frequency, status=status) ! Dont verify since status=-1 is fine here VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%H_order, "H_FILTER_ORDER:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%T_timescale, "T_FILTER_TIMESCALE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%T_frequency, "T_FILTER_FREQUENCY:", rc=status) VERIFY_(status) call Z_Set_Interval(ocean%alarm%Tfilter, ocean%param%filt%T_frequency, status=status) ! Dont verify since status=-1 is fine here VERIFY_(status) call UpdateFromConfig( cf, ocean%param%Filt%T_order, "T_FILTER_ORDER:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%remap_interval, "REMAP_INTERVAL:", rc=status) VERIFY_(status) call Z_Set_Interval(ocean%alarm%remap, ocean%param%remap_interval, status=status) ! Dont verify since status=-1 is fine here VERIFY_(status) call UpdateFromConfig( cf, ocean%param%MassAdvectionOrder, "MASSADVECTIONORDER:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%TracerAdvectionOrder,"TRACERADVECTIONORDER:", rc=status) VERIFY_(status) ! T_Mixed_Layer: call UpdateFromConfig( cf, ocean%ml%do_mixed_layer, "DO_MIXED_LAYER:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%Interval, "ML_INTERVAL:", rc=status) VERIFY_(status) call Z_Set_Interval(ocean%alarm%mixdl, ocean%ml%interval, status=status) ! Dont verify since status=-1 is fine here VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%nLayers, "ML_NLAYERS:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%buffer_layers, "BUFFER_LAYERS:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%dissipat, "ML_DISSIPATION:",rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%penfrac, "ML_PENFRAC:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%pendepth, "ML_PENDEPTH:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%u3coeff, "ML_U3COEFF:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%backU2, "ML_BACKU2:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%sheareff, "ML_SHEAR_EFF:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%ml%qnegdiss, "ML_QNEG_DISS:", rc=status) VERIFY_(status) ! T_HParams: call UpdateFromConfig( cf, ocean%param%hmix%AkH, "AKH:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%AkT, "AKT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%noslip, "USE_NOSLIP:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%laplacian, "USE_LAPLACIAN_VISCOSITY:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%AkU, "AKU:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%Smagorinsky_Lap_Const,"SMAGORINSKY_LAP_CONST:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%Kh_vel_scale,"KH_VEL_SCALE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%biharmonic, "USE_BIHARMONIC_VISCOSITY:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%AaU, "AAU:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%Smagorinsky_Bi_Const,"SMAGORINSKY_BI_CONST:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%Ah_vel_scale,"AH_VEL_SCALE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%hmix%Viscosity_timelimit,"VISCOSITY_TIMELIMIT:", rc=status) VERIFY_(status) ! T_VParams: call UpdateFromConfig( cf, ocean%param%vmix%Interval, "VMIX_INTERVAL:", rc=status) VERIFY_(status) call Z_Set_Interval(ocean%alarm%Vdiff, ocean%param%vmix%interval, status=status) ! Dont verify since status=-1 is fine here VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%CD, "BOTTOM_DRAG_CD:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Extra_bottom_speed, "EXTRA_BOTTOM_SPEED:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Bottom_bl_thickness, "BOTTOM_BL_THICKNESS:", rc=status) VERIFY_(status) ! T_VParams%Kappa_T call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%type_code, "KAPPA_T_TYPE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%max, "KAPPA_T_MAX:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%min, "KAPPA_T_MIN:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%zero_val, "KAPPA_T_0:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%critical_Ri,"KAPPA_T_CRIT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%convect, "KAPPA_T_CONVECT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%power, "KAPPA_T_POWER", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%bott, "KAPPA_T_BOTT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%PGT1, "KAPPA_T_PGT1:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_T%Kraus_coef, "KAPPA_T_KRAUS:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%type_code, "KAPPA_S_TYPE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%max, "KAPPA_S_MAX:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%min, "KAPPA_S_MIN:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%zero_val, "KAPPA_S_0:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%critical_Ri,"KAPPA_S_CRIT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%convect, "KAPPA_S_CONVECT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%power, "KAPPA_S_POWER", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%bott, "KAPPA_S_BOTT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%PGT1, "KAPPA_S_PGT1:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Kappa_S%Kraus_coef, "KAPPA_S_KRAUS:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%type_code, "NU_TYPE:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%max, "NU_MAX:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%min, "NU_MIN:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%zero_val, "NU_0:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%critical_Ri,"NU_CRIT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%convect, "NU_CONVECT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%power, "NU_POWER", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%bott, "NU_BOTT:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%PGT1, "NU_PGT1:", rc=status) VERIFY_(status) call UpdateFromConfig( cf, ocean%param%vmix%Nu%Kraus_coef, "NU_KRAUS:", rc=status) VERIFY_(status) call Update_hmins_from_config(cf,grid,rc=status) VERIFY_(status) binit(:) = ocean%binit(:) call Load_vector_from_config( cf, & LABEL=trim(comp_name)//"_binit:",& VALUES=binit, NFOUND=n, RC=status) if ( n > 0 ) then call write_parallel('Buoyancy targets being changed from restart') ocean%binit(:) = binit(:) else call write_parallel(' no changes to buoyancy targets') endif call write_parallel(' Done with parameters') RETURN_(status) end subroutine Update_Parameters_from_Config subroutine Update_hmins_from_config(cf,grid,rc) type (ESMF_Config) :: cf type (T_Poseidon_Grid) :: grid integer, intent(out), optional :: rc ! Now we need to update hmins from the config file: ! Look for the name in the config file, then loop through the ! values till no more are found. ! Temporaries for reading in and altering hmin, hmax, etc. real :: hmin(grid%kmg), hmax(grid%kmg), hmin_b(grid%kmg) real :: dmin1d(grid%kmg), dmax1d(grid%kmg) real :: default_hmin = 1.0 ! Value to use if hmin <=0 real :: default_hmax = 10000. ! Value to use if hmax <=0 real :: default_hmin_b = 0.001 ! Value to use if hmax <=0 real :: value ! temporary to hold value returned from ConfigGetAttribute integer :: k,km integer :: nhmin,nhmax,nminb,ndmin,ndmax character(ESMF_MAXSTR), parameter :: Iam = 'Update_hmins_from_config' km = grid%kmg hmin(:) = grid%hmin(:) hmax(:) = grid%hmax(:) hmin_b(:) = grid%hmin_b(:) dmin1d(:) = grid%d_min1d(:) dmax1d(:) = grid%d_max1d(:) ! First see if there is a modification to the default values call ESMF_ConfigGetAttribute( cf, VALUE=default_hmin, & LABEL=trim(comp_name)//"_Default_hmin:",& DEFAULT=default_hmin, RC=status) VERIFY_(status) call ESMF_ConfigGetAttribute( cf, VALUE=default_hmax, & LABEL=trim(comp_name)//"_Default_hmax:",& DEFAULT=default_hmax, RC=status) VERIFY_(status) call ESMF_ConfigGetAttribute( cf, VALUE=default_hmin_b, & LABEL=trim(comp_name)//"_Default_hmin_b:",& DEFAULT=default_hmin_b, RC=status) VERIFY_(status) call Load_vector_from_config( cf, & LABEL=trim(comp_name)//"_Hmin:",& VALUES=hmin, NFOUND=nhmin, RC=status) call Load_vector_from_config( cf, & LABEL=trim(comp_name)//"_Hmax:",& VALUES=hmax, NFOUND=nhmax, RC=status) call Load_vector_from_config( cf, & LABEL=trim(comp_name)//"_Hmin_b:",& VALUES=hmin_b, NFOUND=nminb, RC=status) call Load_vector_from_config( cf, & LABEL=trim(comp_name)//"_D_min1d:",& VALUES=dmin1d, NFOUND=ndmin, RC=status) call Load_vector_from_config( cf, & LABEL=trim(comp_name)//"_D_max1d:",& VALUES=dmax1d, NFOUND=ndmax, RC=status) if ( hmin(1) <= 0.0 ) hmin(1) = dmin1d(1) if ( hmin(1) <= 0.0 ) hmin(1) = default_hmin do k=2,km if ( hmin(k) <= 0.0 ) then hmin(k:km) = hmin(k-1) exit endif enddo if ( hmax(1) <= 0.0 ) hmax(1) = dmax1d(1) if ( hmax(1) <= 0.0 ) hmax(1) = default_hmax do k=2,km if ( hmax(k) <= 0.0 ) then hmax(k:km) = hmax(k-1) exit endif enddo if ( hmin_b(1) <= 0.0 ) hmin_b(1) = default_hmin_b do k=2,km if ( hmin_b(k) <= 0.0 ) then hmin_b(k:km) = hmin_b(k-1) exit endif enddo if ( dmin1d(1) <= 0.0 ) dmin1d(1) = hmin(1) do k=1,km if ( dmin1d(k) <= 0.0 ) dmin1d(k) = dmin1d(k-1) + hmin(k) enddo if ( dmax1d(1) <= 0.0 ) dmax1d(1) = hmax(1) do k=1,km if ( dmax1d(k) <= 0.0 ) dmax1d(k) = dmax1d(k-1) + hmax(k) enddo if ( any(grid%hmin(1:km) /= hmin(1:km) ) ) & call write_parallel( 'Hmin is being changed from restart file value') if ( any(grid%hmax(1:km) /= hmax(1:km) ) ) & call write_parallel( 'Hmax is being changed from restart file value') if ( any(grid%hmin_b(1:km) /= hmin_b(1:km) ) ) & call write_parallel( 'Hmin_b is being changed from restart file value') if ( any(grid%d_min1d(1:km) /= dmin1d(1:km) ) ) & call write_parallel( 'D_min1d is being changed from restart file value') if ( any(grid%d_max1d(1:km) /= dmax1d(1:km) ) ) & call write_parallel( 'D_max1d is being changed from restart file value') call write_parallel('Retting grid hmin, hmax values') call Reset_grid_hmins( grid, hmin=hmin,hmax=hmax,hmin_b=hmin_b,& d_min1d=dmin1d,d_max1d=dmax1d) status=ESMF_SUCCESS RETURN_(status) end subroutine Update_hmins_from_config ! This routine looks at a config file subroutine Load_vector_from_config( cf, label, values, nfound, rc ) type (ESMF_Config) :: cf character(*),intent(in) :: label real, intent(inout) :: values(:) integer, optional, intent(out) :: nfound integer, optional, intent(out) :: rc integer :: km integer :: k integer :: kfound real :: value km = size(values) kfound = 0 call ESMF_ConfigFindLabel( cf, LABEL=trim(label), RC=status) if ( present(rc) ) rc=status if ( status == ESMF_SUCCESS ) then do k=1,km call ESMF_ConfigGetAttribute( cf, value, rc=status ) if ( status .ne. ESMF_SUCCESS ) exit kfound = k values(k) = value enddo endif if ( present(nfound) ) nfound = kfound rc = ESMF_SUCCESS end subroutine Load_vector_from_config subroutine Setup_clocks( rc ) integer, intent(out), optional :: rc integer :: reset_ocean_clock type (Z_DateTime) :: ZdateTime type (ESMF_Time) :: Edate type (ESMF_TimeInterval) :: eDt integer :: yy,mm,dd,h,m,s ! If the user wants to reset the ocean clock to be equal to the ! master clock, the config file will have RESET_OCEAN_CLOCK: 1 ! In this case, we reset the time on the ocean clock to be that ! of the master. ! >>>>>>>WARNING: There is a potential mismatch of calendars.<<<<<<<< ! Poseidon uses a gregorian calendar reset_ocean_clock = 0 call UpdateFromConfig( cf, reset_ocean_clock, "RESET_OCEAN_CLOCK:", rc=status) VERIFY_(status) if ( reset_ocean_clock > 0 ) then call write_parallel("Resetting ZEUS clock to match ESMF") call ESMF_ClockGet( clock, currTime = Edate, timeStep = eDt, rc=STATUS) VERIFY_( status ) ! There is a trick that Atanas uses to move the clock back one step ! before initializing, then forward after. We need to add the ! tick here: Edate = Edate + eDt call ESMF_TimeGet( Edate, yy=yy,mm=mm,dd=dd,h=h,m=m,s=s, rc=STATUS) VERIFY_( status ) ZdateTime = Z_as_DateTime(yy,mm,dd,h,m,s) call Z_Set_Time( zclock, ZdateTime, status=status) VERIFY_( status ) call Z_Set_Time( ocean%Alarm%step, ZdateTime , status=status) VERIFY_( status ) else call write_parallel("Not altering zeus clock time") ZdateTime = Z_as_DateTime( zclock ) endif ! create an alarm for stopping each step. Give it a ring interval ! of 1.0 so that it will ring every time if we fail to set it ! properly call Z_Init_Alarm( zclock, POS_MAPL_internal_state%EOR_Alarm, & & NAME='Poseidon Stop Alarm', & & ORIGIN=ZdateTime, & & INTERVAL=1.0, & & RING=.FALSE. ) RETURN_( status ) end subroutine Setup_clocks subroutine choose_diagnostics(rc) integer, intent(out),optional:: rc character(ESMF_MAXSTR) :: name integer :: len integer :: i call write_parallel(' ') call write_parallel(' SETTING UP OCEAN HISTORY DIAGNOSTICS') ! Get the diagnostics from the config file ! If no History Diagnostics: is in the config file, skip this ! without causing a problem. len = ESMF_ConfigGetLen( cf, trim(comp_name)//"_HISTORY_DIAGNOSTICS:", rc = status ) if( status == ESMF_SUCCESS ) then call ESMF_ConfigFindLabel( cf, trim(comp_name)//"_HISTORY_DIAGNOSTICS:", rc=status) VERIFY_(status) do i=1,len call ESMF_ConfigGetAttribute( cf, name, rc=status ) if ( status .ne. ESMF_SUCCESS ) exit call MarkDiagForHistory( ocean%diag, name=trim(name) ) call write_parallel(' Adding Poseidon history diagnostic '//trim(name) ) enddo else call write_parallel('"'//trim(comp_name)//'_HISTORY_DIAGNOSTICS:" not found in configuration') call write_parallel(' default history quantities will be used') endif call PrintWanted( ocean%diag ) if ( present(rc) ) rc=ESMF_SUCCESS end subroutine choose_diagnostics subroutine Setup_history( rc ) integer, intent(out),optional:: rc integer :: input(2) integer :: years integer :: months integer :: days integer :: hours integer :: minutes integer :: seconds type (ESMF_TimeInterval) :: dt type (Z_DateTime) :: dateTime, dateTime2 type (Z_TimeInterval) :: timeInterval character(ESMF_MAXSTR) :: filename,template,ctl_name status = 0 POS_MAPL_internal_state%write_control_file = .true. ! history and history accumulator alarms are handled at the ! POS_Plug level - not witin poseidon call ESMF_ConfigGetAttribute( cf, years, & LABEL=trim(comp_name)//"_HISTORY_YEARS:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, months, & LABEL=trim(comp_name)//"_HISTORY_MONTHS:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, days, & LABEL=trim(comp_name)//"_HISTORY_DAYS:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, hours, & LABEL=trim(comp_name)//"_HISTORY_HOURS:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, minutes, & LABEL=trim(comp_name)//"_HISTORY_MINUTES:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, seconds, & LABEL=trim(comp_name)//"_HISTORY_SECONDS:", & DEFAULT=0, RC=status) VERIFY_(status) dateTime%year=1 dateTime%month=1 dateTime%day=1 dateTime%hour=0 dateTime%minute=0 dateTime%second=0 if ( months > 0 .or. years > 0 ) then ! Set the alarm origin to the present time plus N months or years dateTime = Z_as_DateTime( zclock ) dateTime%month = dateTime%month + months dateTime%Year = dateTime%year + years timeInterval = Z_as_TimeInterval( SECONDS=1.) else ! for all others, the origin is 01/01/0001 00:00:00 dateTime2=dateTime dateTime2%second = dateTime2%second + seconds dateTime2%minute = dateTime2%minute + minutes dateTime2%hour = dateTime2%hour + hours dateTime2%day = dateTime2%day + days timeInterval = dateTime2 - DateTime endif call Z_Init_Alarm( zclock, POS_MAPL_internal_state%history_Alarm, & & NAME='Poseidon History Alarm', & & ORIGIN=dateTime, & & INTERVAL=timeInterval, & & RING=.FALSE. ) ! History is sampled more frequently than it is written. ! Here we set the sampling interval call ESMF_ConfigGetAttribute( cf, years, & LABEL=trim(comp_name)//"_HISTORY_SAMPLING_YEARS:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, months, & LABEL=trim(comp_name)//"_HISTORY_SAMPLING_MONTHS:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, days, & LABEL=trim(comp_name)//"_HISTORY_SAMPLING_DAYS:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, hours, & LABEL=trim(comp_name)//"_HISTORY_SAMPLING_HOURS:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, minutes, & LABEL=trim(comp_name)//"_HISTORY_SAMPLING_MINUTES:", & DEFAULT=0, RC=status) call ESMF_ConfigGetAttribute( cf, seconds, & LABEL=trim(comp_name)//"_HISTORY_SAMPLING_SECONDS:", & DEFAULT=0, RC=status) VERIFY_(status) dateTime%year=1 dateTime%month=1 dateTime%day=1 dateTime%hour=0 dateTime%minute=0 dateTime%second=0 if ( months > 0 .or. years > 0 ) then ! Set the alarm origin to the present time plus N months or years dateTime = Z_as_DateTime( zclock ) dateTime%month = dateTime%month + months dateTime%Year = dateTime%year + years timeInterval = Z_as_TimeInterval( SECONDS=1.) else ! for all others, the origin is 01/01/0001 00:00:00 dateTime2=dateTime dateTime2%second = dateTime2%second + seconds dateTime2%minute = dateTime2%minute + minutes dateTime2%hour = dateTime2%hour + hours dateTime2%day = dateTime2%day + days timeInterval = dateTime2 - DateTime endif call Z_Init_Alarm( zclock, POS_MAPL_internal_state%history_Acc_Alarm, & & NAME='Poseidon History Sampling Alarm', & & ORIGIN=dateTime, & & INTERVAL=timeInterval, & & RING=.FALSE. ) ! It is expected that diagnostics have already been ! chosen. call write_parallel(' ') call write_parallel(' CREATING AN OCEAN HISTORY OBJECT') call Create_ohist( ocean, oforce, ohist ) call ESMF_ConfigGetAttribute( cf, filename, & LABEL=trim(comp_name)//"_HISTORY_FILENAME:", & DEFAULT='o_hist', RC=status) VERIFY_(status) POS_MAPL_internal_state%history_unit = Getfile( trim(filename),& FORM='UNFORMATTED',RC=status) VERIFY_(status) ! Use the history filename as the template, but pre-pend '^' ! as a default. call ESMF_ConfigGetAttribute( cf, template, & LABEL=trim(comp_name)//"_HISTORY_GRADS_TEMPLATE:", & DEFAULT='^'//trim(filename), RC=status ) VERIFY_(status) POS_MAPL_internal_state%ohist_template = trim(template) call ESMF_ConfigGetAttribute( cf, ctl_name, & LABEL=trim(comp_name)//"_HISTORY_GRADS_FILENAME:", & DEFAULT=trim(filename)//'.ctl', RC=status ) VERIFY_(status) POS_MAPL_internal_state%grads_descriptor_unit = Getfile( trim(ctl_name), & FORM='FORMATTED',RC=status) VERIFY_(status) call write_parallel(' Ocean history file: '//trim(filename) ) call write_parallel(' Grads control file: '//trim(ctl_name) ) call write_parallel(' Grads file template: '//trim(template) ) call write_parallel(' Ocean History creation completed') RETURN_(status) end subroutine Setup_history end subroutine Initialize !================================================================================= !BOP ! !IROUTINE: Run -- Run method for POS_Plug wrapper ! !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 supervisor clock integer, optional, intent( OUT) :: rc ! Error code: ! !DESCRIPTION: ! This routine takes the following steps: ! \begin{enumerate} ! \item Extract the private Poseidon composite state from the gridded component, ! and associate poseidon style pointers with the internals. ! \item Extract the import and export pointers, then copy forcing from ! the import state to the Poseidon forcing object. Do units conversion ! and vector rotations here, as well as compute some composite ! quantities (buoyancy flux) from the defined imports. ! \item Set up the private Zeus clock to enable Poseidon to run until the end ! of the controlling time step. This is simply done by making a ``stop\_alarm'' ! which will ring at the current (ESMF) time plus one (ESMF) time step. ! \item Update the ocean state with the penetrating radiative heating ! \item Step Poseidon through its timesteps until ``stop\_alarm'' rings. ! \item During this stepping, accumulate private Poseidon history and ! write it out, if required. Use of Poseidon private history ! is retained for backward compatibility, but use of MAPL history is ! now preferred. (The two are NOT mutually exclusive, and both work). ! \item After running for the required time, update the export state quantities ! \item Exit ! \end{enumerate} !EOP ! Locals integer :: IM, JM, LM integer :: IMw, JMw integer :: num_ocean_calls = 1 logical :: ocean_seg_start = .false. logical :: ocean_seg_end = .false. real, pointer :: TAUX(:,:) real, pointer :: TAUY(:,:) real, pointer :: U3(:,:) real, pointer :: PS (:,:) real, pointer :: HEAT(:,:,:) real, pointer :: HeatFlux(:,:) real, pointer :: WaterFlux(:,:) real, pointer :: SaltFlux(:,:) real, pointer :: RadHeating(:,:,:) real, pointer :: TS (:,:) real, pointer :: SS (:,:) real, pointer :: US (:,:) real, pointer :: VS (:,:) real, pointer :: DH (:,:,:) real, pointer :: MASK(:,:,:) type(MAPL_MetaComp), pointer :: MAPL type(POS_MAPL_Type), pointer :: POS_MAPL_internal_state type(POS_MAPLWrap_Type) :: wrap type (T_poseidon_ocean ), pointer :: ocean type (T_poseidon_grid ), pointer :: grid type (T_poseidon_forcing), pointer :: oforce type (T_poseidon_history), pointer :: ohist type (Z_clock), pointer :: zclock type (Z_alarm), pointer :: stop_Alarm type (Z_alarm), pointer :: history_Acc_Alarm type (Z_alarm), pointer :: history_Alarm type(ESMF_Time) :: eCurrTime type(ESMF_TimeInterval) :: eTimeStep type(ESMF_Time) :: eEndTime type(Z_DateTime) :: zCurrTime type(Z_DateTime) :: zEndTime type(Z_TimeInterval) :: zTimeStep type(Z_DateTime) :: histDateTime type(Z_DateTimeInterval) :: dti integer :: i1,i2,j1,j2,km ! Begin !------ ! Get the component's name and set-up traceback handle. ! ----------------------------------------------------- Iam = "Run" call ESMF_GridCompGet( gc, NAME=comp_name, RC=status ) VERIFY_(status) Iam = trim(comp_name) // Iam ! Get my internal MAPL_Generic state !----------------------------------- call MAPL_GetObjectFromGC ( GC, MAPL, RC=status) VERIFY_(STATUS) ! Profilers !---------- call MAPL_TimerOn (MAPL,"TOTAL") call MAPL_TimerOn (MAPL,"RUN" ) ! Get GuestModel's private internal state !--------------------------------- CALL ESMF_UserCompGetInternalState( GC, 'POS_MAPL_state', WRAP, STATUS ) VERIFY_(STATUS) POS_MAPL_internal_state => WRAP%PTR ! shortcut pointers !--------------------------------- ocean => POS_MAPL_internal_state%ocean grid => POS_MAPL_internal_state%grid oforce => POS_MAPL_internal_state%oforce ohist => POS_MAPL_internal_state%ohist zclock => POS_MAPL_internal_state%zclock stop_Alarm => POS_MAPL_internal_state%EOR_Alarm history_Alarm => POS_MAPL_internal_state%history_Alarm history_Acc_Alarm => POS_MAPL_internal_state%history_Acc_Alarm i1 = grid%i1 i2 = grid%i2 j1 = grid%j1 j2 = grid%j2 km = grid%kmg ! Get IMPORT pointers !-------------------- call MAPL_GetPointer(IMPORT, TAUX, 'TAUX', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT, TAUY, 'TAUY', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT, U3, 'USTAR3', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT, PS, 'PS' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT, HeatFlux, 'HFLX', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT, WaterFlux, 'QFLX', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT, SaltFlux, 'SFLX', RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(IMPORT, RadHeating, 'SWHEAT', RC=STATUS); VERIFY_(STATUS) ! Get EXPORT pointers !-------------------- call MAPL_GetPointer(EXPORT, US, 'US' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, VS, 'VS' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, TS, 'TS' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, SS, 'SS' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, DH, 'DH' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(EXPORT, MASK, 'MASK' , RC=STATUS); VERIFY_(STATUS) ! We are ignoring surface pressure for now PSS 3/9/07 ! this is because grottoes come across as PS = 1.E+15 !call assign_forcing( grid, PS, oforce%P0, rc=STATUS) !VERIFY_(STATUS) oforce%P0 = 100000.0 call write_parallel(' Starting to assign forcing ') call assign_forcing( grid, HeatFlux, oforce%HeatFlux, rc=STATUS) VERIFY_(STATUS) call assign_forcing( grid, WaterFlux, oforce%WaterFlux, rc=STATUS) VERIFY_(STATUS) ! We get surface heat flux with the enthalpy of rain water ! computed as waterFlux * T * C_p, where T is in Kelvin ! but we use it as though T is in Centigrade ! remove the difference before Poseidon uses it: oforce%HeatFlux = oforce%HeatFlux - oforce%WaterFlux * CTOKELV * CP_OCN call assign_forcing( grid, SaltFlux, oforce%SaltFlux, rc=STATUS) VERIFY_(STATUS) call assign_vector_forcing( grid, TAUX, TAUY, oforce%taux, oforce%tauy, rc=STATUS) VERIFY_(STATUS) call assign_forcing( grid, U3, oforce%U3, rc=STATUS) VERIFY_(STATUS) call Compute_Buoyancy_Flux( grid, oforce%WaterFlux, & oforce%HeatFlux,oforce%SaltFlux,& ocean%state%T(:,:,1),& ocean%state%S(:,:,1)+ocean%state%Sbias, & ocean%state%H(:,:,1), & oforce%BuoyancyFlux, rc=status) VERIFY_(STATUS) oforce%Qrad_net(:,:) = 0.0 ! We handle this directly, not in poseidon ! Run Poseidon !-------------------------- ! get the ESMF clock time and DT ! get the zeus clock time and dt ! run poseidon til you get up to or past call ESMF_ClockGet( clock, currTime=eCurrTime, & timeStep=eTimeStep, rc=STATUS) VERIFY_(STATUS) eEndTime = eCurrTime + eTimeStep zEndTime = eEndTime zCurrTime = Z_as_DateTime( zclock ) ! Set the stop_alarm to start ringing when it reaches ! the end time. (Set its origin to zEndTime, and its ! interval to 1.0) call Z_Set_Time( stop_alarm, zEndTime , status=STATUS) VERIFY_(STATUS) call Z_Set_Interval( stop_alarm, 1.0 , status=STATUS) VERIFY_(STATUS) call Z_Shut_Alarm( stop_alarm, status=STATUS) VERIFY_(STATUS) call Apply_radiative_heating( eTimeStep, rc=STATUS) do while ( .not. Z_Ringing( stop_Alarm ) ) ! Poseidon uses a tick-then-run paradigm call Z_tick( zclock ) if ( Z_Ringing( history_Acc_Alarm ) .or. & Z_Ringing( history_Alarm ) ) then call EnableDiagnostics( ocean%diag) else call DisableDiagnostics( ocean%diag) endif call MAPL_TimerOn (MAPL,"--PosRun" ) call run_poseidon( ocean,oforce ) call MAPL_TimerOff(MAPL,"--PosRun" ) ! if ( use_sponge .and. Z_RINGING( sponge_alarm ) ) call apply_sponge if ( Z_RINGING( history_Acc_Alarm) ) then call MAPL_TimerOn (MAPL,"--PosHistory" ) call Accumulate_Ohist( ohist, 1.0 ) call Z_Shut_alarm( history_Acc_Alarm ) call ResetDiagnostic( ocean%diag ) call MAPL_TimerOff(MAPL,"--PosHistory" ) endif if ( Z_Ringing( history_Alarm ) ) then call MAPL_TimerOn (MAPL,"--PosHistory" ) call write_parallel(' WRITING OCEAN HISTORY') call Z_Shut_Alarm( history_Alarm ) dti = POS_MAPL_internal_state%historyInterval call Write_ohist( ohist, & POS_MAPL_internal_state%history_unit, & POS_MAPL_internal_state%grads_descriptor_unit , & POS_MAPL_internal_state%write_control_file,& n_secs = int(dti%seconds) ,& n_mins = int(dti%minutes) ,& n_hours= int(dti%hours) ,& n_days = int(dti%days) ,& n_months=int(dti%months) ,& n_years= int(dti%years) ,& file_template=TRIM(POS_MAPL_internal_state%ohist_template)) POS_MAPL_internal_state%write_control_file = .false. call ResetDiagnostic( ocean%diag ) call Initialize_ohist(ohist) ! If the history alarm increment is not a multiple ! of seconds (or hrs,mins, days), we need to reset ! origin for months and years if ( dti%months > 0 .or. dti%years > 0 ) then histDateTime = Z_as_DateTime( history_Alarm ) histDateTime%Month = histDateTime%Month + dti%months histDateTime%Year = histDateTime%Year + dti%years call Z_Set_Time(history_Alarm,histDateTime) endif call MAPL_TimerOff(MAPL,"--PosHistory" ) endif zCurrTime = Z_as_DateTime( zclock ) enddo ! Get the A grid currents if(associated(US)) US = ocean%state%U(i1:i2,j1:j2,1)*grid%cos_chi(i1:i2,j1:j2) & + ocean%state%V(i1:i2,j1:j2,1)*grid%sin_chi(i1:i2,j1:j2) if(associated(VS)) VS = ocean%state%V(i1:i2,j1:j2,1)*grid%cos_chi(i1:i2,j1:j2) & - ocean%state%U(i1:i2,j1:j2,1)*grid%sin_chi(i1:i2,j1:j2) if(associated(TS )) then TS = ocean%state%T(i1:i2,j1:j2,1) + CTOKELV end if if(associated(SS )) then SS = ocean%state%S(i1:i2,j1:j2,1) + ocean%state%sbias end if if(associated(DH)) then DH = ocean%state%H(i1:i2,j1:j2,:) end if if(associated(MASK)) then MASK = grid%bh3(i1:i2,j1:j2,:) end if call MAPL_TimerOff(MAPL,"RUN" ) call MAPL_TimerOff(MAPL,"TOTAL" ) ! All Done !--------- RETURN_(ESMF_SUCCESS) contains ! Take data from an import state (such as fluxes) ! and put it on a local copy that is associated ! with the poseidon grid. ! Multiply by a conversion factor if needed ! Halo fill the result subroutine assign_forcing( grid, input, output, conversion, rc) type (T_poseidon_grid), intent(in) :: grid real, intent(in) :: input(:,:) real, intent(out) :: output(:,:) real, intent(in), optional :: conversion integer, intent(out), optional :: rc real :: mult logical :: ok integer :: i1,i2,j1,j2 call PGRID_findcore( grid, output, i1,i2,j1,j2, ok ) ASSERT_( ok ) output( i1:i2, j1:j2 ) = input(:,:) if ( present(conversion) ) then output( i1:i2,j1:j2) = output( i1:i2,j1:j2 ) * conversion endif ! Fill halos call HaloFill(grid%news, output, tripole_stagger=gmA_STAGGER, vector=.false.) if ( present(rc) ) rc = ESMF_SUCCESS end subroutine assign_forcing subroutine assign_vector_forcing( grid, in_E, in_N, out_u, out_v, conversion, rc) type (T_poseidon_grid), intent(in) :: grid real, intent(in) :: in_E(:,:) real, intent(in) :: in_N(:,:) real, intent(out) :: out_u(:,:) real, intent(out) :: out_v(:,:) real, intent(in), optional :: conversion integer, intent(out), optional :: rc real :: mult logical :: ok integer :: i1,i2,j1,j2 ! If you want to print out a specific spot on the world grid ! you set itest and jtest according to world coordinates. ! then if (do_report), the present PE has the right value integer :: itest,jtest logical :: do_report !------------------------------------------------------------------ ! If you want to print out a specific spot on the world grid ! you set itest and jtest according to world coordinates. ! then if (do_report), the present PE has the right value itest=203-grid%i_world_offset jtest=150-grid%j_world_offset do_report = ( itest > 2 .and. itest < grid%img-2 .and. & jtest > 2 .and. jtest < grid%jmg-2 ) call PGRID_findcore( grid, out_u, i1,i2,j1,j2, ok ) ASSERT_( ok ) out_u( i1:i2, j1:j2 ) = in_E(:,:)*grid%cos_chi(i1:i2,j1:j2) & - in_N(:,:)*grid%sin_chi(i1:i2,j1:j2) out_v( i1:i2, j1:j2 ) = in_N(:,:)*grid%cos_chi(i1:i2,j1:j2) & + in_E(:,:)*grid%sin_chi(i1:i2,j1:j2) if ( present(conversion) ) then out_u( i1:i2,j1:j2) = out_u( i1:i2,j1:j2 ) * conversion out_v( i1:i2,j1:j2) = out_v( i1:i2,j1:j2 ) * conversion endif ! Fill halos call HaloFill( grid%news, out_u, tripole_stagger=gmA_STAGGER, vector=.true. ) call HaloFill( grid%news, out_v, tripole_stagger=gmA_STAGGER, vector=.true. ) if ( present(rc) ) rc = ESMF_SUCCESS end subroutine assign_vector_forcing subroutine Compute_Buoyancy_Flux( g, Waterflux, Heatflux, Saltflux,& SST, SSS, H, BuoyFlux, rc ) type (T_poseidon_grid), intent(in) :: g real, intent(in) :: Waterflux(:,:) ! Water flux (kg/m^2/s) into ocean real, intent(in) :: Heatflux(:,:) ! Enthalpy Flux (W/m^2) into ocean real, intent(in) :: Saltflux(:,:) ! Salt Flux (kg/m^2/s) into ocean real, intent(in) :: SST(:,:) ! SST (C) real, intent(in) :: SSS(:,:) ! SSS (psu) real, intent(in) :: H(:,:) ! mixed layer mass (dyn-m) real, intent(out) :: BuoyFlux(:,:) integer, intent(out),optional :: rc real :: dTdt(g%img,g%jmg),dSdt(g%img,g%jmg),Tflux(g%img,g%jmg) real :: Hr(g%img,g%jmg) character(ESMF_MAXSTR) :: message integer :: i,j integer :: itest,jtest logical :: do_report itest=1-g%i_world_offset jtest=50-g%j_world_offset do_report = ( itest > 2 .and. itest < g%img-2 .and. & jtest > 2 .and. jtest < g%jmg-2 ) ! Compute the rate of change in temperature ! This includes the addition of mass with some temperature ! (bundled in Heatflux) Tflux(:,:) = Heatflux(:,:) * RHOCP_R - Waterflux(:,:) * SST(:,:) * RHO_OCN_R where ( H(:,:) > 0.0 ) Hr(:,:) = 1./H(:,:) elsewhere Hr(:,:) = 0.0 endwhere dTdt(:,:) = Tflux(:,:) * Hr(:,:) dSdt(:,:) = ( Saltflux(:,:) * 1000. & - Waterflux(:,:) * RHO_OCN_R * SSS(:,:) ) * Hr(:,:) do j=1,g%jmg do i=1,g%img if( Hr(i,j) > 0.0 ) then BuoyFlux(i,j) = ( dbdtemp(SST(i,j),SSS(i,j))*dTdt(i,j) & + dbdsalt(SST(i,j),SSS(i,j))*dSdt(i,j) ) * H(i,j) else BuoyFlux(i,j) = 0.0 endif enddo enddo if ( present(rc) ) rc = ESMF_SUCCESS end subroutine Compute_Buoyancy_Flux subroutine Apply_radiative_heating( dt, RC ) type (ESMF_TimeInterval), intent(in) :: dt integer, intent(out),optional :: RC integer :: idt real :: delt call ESMF_TimeIntervalGet( dt,s=idt,rc=STATUS) VERIFY_(status) delt = idt ASSERT_( SIZE(RadHeating,1)==(i2-i1+1) ) ASSERT_( SIZE(RadHeating,2)==(j2-j1+1) ) ASSERT_( SIZE(RadHeating,3)==km ) ! tell the model that we have done radiation already oforce%use_external_radiation_routine = .true. ! Assign the RadHeating field in oforce (for diagnostic only) oforce%RadHeating(i1:i2,j1:j2,1:km) = RadHeating(:,:,:) where ( ocean%state%H(i1:i2,j1:j2,1:km) > 0.01 ) & ocean%state%T(i1:i2,j1:j2,1:km) = & ocean%state%T(i1:i2,j1:j2,1:km) & + RadHeating(:,:,:) * RHOCP_R * delt / & ocean%state%H(i1:i2,j1:j2,1:km) ! Fill halos call HaloFill( grid%news, ocean%state%T, tripole_stagger=gmA_STAGGER, vector=.false. ) call HaloFill( grid%news, oforce%RadHeating, tripole_stagger=gmA_STAGGER, vector=.false. ) end subroutine Apply_radiative_heating end subroutine Run !BOP ! !IROUTINE: Finalize -- Finalize method for POS_Plug wrapper ! !INTERFACE: subroutine Finalize ( 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 supervisor clock integer, optional, intent( OUT) :: rc ! Error code: ! !DESCRIPTION: ! Finalize performs the following tasks: ! ! \begin{enumerate} ! \item Extract the private Poseidon composite state from the component ! \item Determine the requested name and type of the private Poseidon checkpoint file ! \item Save the checkpoint data. At this time, only the private ZDF format is ! allowed. ! \item Deallocate any space used for history, clocks and the poseidon state. ! \item Call MAPL\_GenericFinalize for full clean-up ! \end{enumerate} !EOP type(MAPL_MetaComp), pointer :: State type(ESMF_Config) :: CF type(POS_MAPL_Type), pointer :: POS_MAPL_internal_state type(POS_MAPLWrap_Type) :: wrap type(T_Poseidon_ocean), pointer :: ocean type(T_Poseidon_grid), pointer :: grid integer :: iu ! file unit number character(ESMF_MAXSTR) :: filename character(ESMF_MAXSTR) :: filetype ! Get the target components name and set-up traceback handle. ! ----------------------------------------------------------- Iam = "Finalize" call ESMF_GridCompGet( gc, NAME=comp_name, RC=status ) VERIFY_(STATUS) Iam = trim(comp_name) // Iam ! Get my internal MAPL_Generic state !----------------------------------- call MAPL_GetObjectFromGC ( GC, STATE, RC=status) VERIFY_(STATUS) ! Profilers !---------- call MAPL_TimerOn(STATE,"TOTAL" ) call MAPL_TimerOn(STATE,"FINALIZE") ! Get the grid !---------------------------- call ESMF_GridCompGet( GC, CONFIG = CF, RC=status ) VERIFY_(STATUS) ! Get the Plug's private internal state !-------------------------------------- CALL ESMF_UserCompGetInternalState( GC, 'POS_MAPL_state', WRAP, STATUS ) VERIFY_(STATUS) POS_MAPL_internal_state => wrap%ptr ocean => POS_MAPL_internal_state%ocean ! Get the name of the Poseidon old style restart file call ESMF_ConfigGetAttribute( cf, VALUE=filetype, & LABEL=trim(comp_name)//"_POSEIDON_CHECKPOINT_TYPE:", & DEFAULT='zdf',RC=status) VERIFY_(status) call ESMF_ConfigGetAttribute( cf, filename, & LABEL=trim(comp_name)//"_POSEIDON_CHECKPOINT_FILE:", & DEFAULT='none',RC=status) VERIFY_(status) if ( filetype == 'zdf' .or. filetype == 'ZDF' ) then if ( filename .ne. 'none' ) then iu = Getfile( trim(filename), form='unformatted', RC=status) VERIFY_(status) call Save_restart( iu, ocean ) call Free_File(iu) else call write_parallel('No POSEIDON_CHECKPOINT_FILE: found in config file') call write_parallel(' not writing checkpoint') endif else call write_parallel('Only ZDF file type is supported for Poseidon as of now') call write_parallel(' Going to write zdf anyway....') !!!!!!!!!!!!!!!!!! Temporary patch to force zdf restart if ( filename .ne. 'none' ) then iu = Getfile( trim(filename), form='unformatted', RC=status) VERIFY_(status) call Save_restart( iu, ocean ) call Free_File(iu) else call write_parallel('No POSEIDON_CHECKPOINT_FILE: found in config file') call write_parallel(' not writing checkpoint') endif endif call Finalize_history( rc=STATUS ) VERIFY_(STATUS) call Finalize_clocks( rc=STATUS ) VERIFY_(STATUS) call Destroy_Poseidon( POS_MAPL_internal_state , rc=STATUS) VERIFY_(STATUS) !!!!!!!!!!!!!!!!!! call MAPL_TimerOff(STATE,"FINALIZE") call MAPL_TimerOff(STATE,"TOTAL" ) ! Generic initialize ! ------------------ call MAPL_GenericFinalize( GC, IMPORT, EXPORT, CLOCK, RC=status ) VERIFY_(STATUS) ! All Done !--------- RETURN_(ESMF_SUCCESS) contains subroutine Destroy_Poseidon(p,rc) type(POS_MAPL_Type) :: p integer, optional, intent(out) :: rc call Destroy_ocean( p%ocean ) call Destroy_forcing( p%oforce ) call Destroy_ohist( p%ohist ) call Destroy_grid( p%grid ) call Z_Free_clock( p%zclock ) RETURN_(ESMF_SUCCESS) end subroutine Destroy_Poseidon subroutine Finalize_history( rc ) integer, intent(out), optional :: rc call Z_Free_Alarm( POS_MAPL_internal_state%history_Alarm ) call Z_Free_Alarm( POS_MAPL_internal_state%history_Acc_Alarm ) RETURN_(status) end subroutine Finalize_history subroutine Finalize_clocks( rc ) integer, intent(out), optional :: rc call Z_Free_Alarm( POS_MAPL_internal_state%EOR_Alarm ) call Z_Free_Clock( POS_MAPL_internal_state%zclock ) RETURN_(status) end subroutine Finalize_clocks end subroutine Finalize ! Utility routines subroutine UpdateRealFromConfig( cf, val, name, rc ) type (ESMF_Config) :: cf real, intent(inout) :: val character(*), intent(in) :: name integer,optional,intent(out) :: rc real :: tmp integer :: status call ESMF_ConfigGetAttribute( cf, VALUE=tmp, & LABEL=trim(comp_name)//"_"//trim(name),& DEFAULT=val, RC=status) VERIFY_(status) if ( val /= tmp ) then call WRITE_PARALLEL( 'Changing Poseidon Parameter "' // & Trim(name)//'" via Config file') val = tmp endif RETURN_(status) end subroutine UpdateRealFromConfig subroutine UpdateIntFromConfig( cf, val, name, rc ) type (ESMF_Config) :: cf integer, intent(inout) :: val character(*), intent(in) :: name integer,optional,intent(out) :: rc integer :: tmp integer :: status call ESMF_ConfigGetAttribute( cf, VALUE=tmp, & LABEL=trim(comp_name)//"_"//trim(name),& DEFAULT=val, RC=status) VERIFY_(status) if ( val /= tmp ) then call WRITE_PARALLEL( 'Changing Poseidon Parameter "' // & Trim(name)//'" via Config file') val = tmp endif RETURN_(status) end subroutine UpdateIntFromConfig subroutine UpdateLogFromConfig( cf, val, name, rc ) type (ESMF_Config) :: cf logical, intent(inout) :: val character(*), intent(in) :: name integer,optional,intent(out) :: rc integer :: tmp,tmpval integer :: status tmpval = 0 if (val) tmpval=1 call ESMF_ConfigGetAttribute( cf, VALUE=tmp, & LABEL=trim(comp_name)//"_"//trim(name),& DEFAULT=tmpval, RC=status) VERIFY_(status) if ( tmpval /= tmp ) then call WRITE_PARALLEL( 'Changing Poseidon Parameter "' // & Trim(name)//'" via Config file') val = ( tmp /= 0 ) endif RETURN_(status) end subroutine UpdateLogFromConfig ! Convert a Zeus time interval into an ESMF time interval ! Note: zeus time intervals are stored in default real seconds, ! but ESMF time intervals do not accept real seconds. ! This is written as an elemental subroutine so that ! it can be used in assignment operations. subroutine Z_to_E_TimeInterval( Einterval, Zinterval) type (ESMF_TimeInterval), intent(inout) :: Einterval ! ESMF TimeInterval out type (Z_TimeInterval), intent(in) :: Zinterval ! Zeus TimeInterval in integer :: nseconds nseconds = Z_as_Seconds(Zinterval) call ESMF_TimeIntervalSet( Einterval, s=nseconds ) end subroutine Z_to_E_TimeInterval ! Convert an ESMF time interval into a Zeus time interval ! Note: zeus time intervals are stored in default real seconds, ! but ESMF time intervals do not accept real seconds. ! This is written as an elemental subroutine so that ! it can be used in assignment operations. subroutine E_to_Z_TimeInterval( Zinterval, Einterval) type (Z_TimeInterval), intent(inout) :: Zinterval ! Zeus TimeInterval out type (ESMF_TimeInterval), intent(in) :: Einterval ! ESMF TimeInterval in integer :: nseconds call ESMF_TimeIntervalGet( Einterval, s=nseconds ) Zinterval = Z_as_TimeInterval(seconds=nseconds) end subroutine E_to_Z_TimeInterval ! Convert an ESMF Time into a Zeus DateTime subroutine E_to_Z_DateTime( Zdate, Edate) type (Z_DateTime), intent(inout) :: Zdate ! Zeus dateTime out type (ESMF_Time), intent(in) :: Edate ! ESMF time in integer :: yy,mm,dd,h,m,s type (ESMF_calendarType) :: ctype call ESMF_TimeGet( Edate, yy=yy,mm=mm,dd=dd,h=h,m=m,s=s,calendarType=ctype) Zdate = Z_as_DateTime(yy,mm,dd,h,m,s) end subroutine E_to_Z_DateTime ! Convert an ESMF Time into a Zeus DateTime subroutine Z_to_E_DateTime( Edate, Zdate) type (ESMF_Time), intent(inout) :: Edate ! ESMF time out type (Z_DateTime), intent(in) :: Zdate ! Zeus dateTime in call ESMF_TimeSet( Edate, yy=Zdate%year,mm=Zdate%month,dd=Zdate%day,& h=Zdate%hour,m=Zdate%minute,s=int(Zdate%second) ) end subroutine Z_to_E_DateTime !==================================================================== end module POS_GEOS5PlugMod