C $Header: /u/gcmpack/MITgcm/model/src/ini_grid.F,v 1.35 2013/02/17 02:18:16 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP C !ROUTINE: INI_GRID C !INTERFACE: SUBROUTINE INI_GRID( myThid ) C !DESCRIPTION: C These arrays are used throughout the code in evaluating gradients, C integrals and spatial avarages. This routine is called separately C by each thread and initializes only the region of the domain it is C "responsible" for. C !CALLING SEQUENCE: C INI_GRID C | -- LOAD_GRID_SPACING C | -- INI_VERTICAL_GRID C | / INI_CARTESIAN_GRID C | / INI_SPHERICAL_POLAR_GRID C | \ INI_CURVILINEAR_GRID C | \ INI_CYLINDER_GRID C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #ifdef ALLOW_MNC #include "MNC_PARAMS.h" #endif #ifdef ALLOW_MONITOR #include "MONITOR.h" #endif C !INPUT/OUTPUT PARAMETERS: C myThid :: my Thread Id number INTEGER myThid C !FUNCTIONS: LOGICAL MASTER_CPU_IO EXTERNAL MASTER_CPU_IO C !LOCAL VARIABLES: C bi, bj :: tile indices C i, j :: Loop counters C msgBuf :: Informational/error message buffer INTEGER bi, bj INTEGER i, j CHARACTER*(MAX_LEN_MBUF) msgBuf CEOP C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Load grid spacing (vector) from files CALL LOAD_GRID_SPACING( myThid ) C-- Set up vertical grid and coordinate system CALL INI_VERTICAL_GRID( myThid ) C-- Initialise (everywhere) all horizontal grid array to null value C Note: some arrays are not defined in some parts of the halo C region. We set them to zero here for safety. If they are ever C referred to, especially in the denominator then it is a mistake! DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx xC(i,j,bi,bj) = 0. yC(i,j,bi,bj) = 0. xG(i,j,bi,bj) = 0. yG(i,j,bi,bj) = 0. dxC(i,j,bi,bj) = 0. dyC(i,j,bi,bj) = 0. dxG(i,j,bi,bj) = 0. dyG(i,j,bi,bj) = 0. dxF(i,j,bi,bj) = 0. dyF(i,j,bi,bj) = 0. dxV(i,j,bi,bj) = 0. dyU(i,j,bi,bj) = 0. rA(i,j,bi,bj) = 0. rAz(i,j,bi,bj) = 0. rAw(i,j,bi,bj) = 0. rAs(i,j,bi,bj) = 0. recip_dxG(i,j,bi,bj) = 0. recip_dyG(i,j,bi,bj) = 0. recip_dxC(i,j,bi,bj) = 0. recip_dyC(i,j,bi,bj) = 0. recip_dxF(i,j,bi,bj) = 0. recip_dyF(i,j,bi,bj) = 0. recip_dxV(i,j,bi,bj) = 0. recip_dyU(i,j,bi,bj) = 0. recip_rA (i,j,bi,bj) = 0. recip_rAs(i,j,bi,bj) = 0. recip_rAw(i,j,bi,bj) = 0. recip_rAz(i,j,bi,bj) = 0. tanPhiAtU(i,j,bi,bj) = 0. tanPhiAtV(i,j,bi,bj) = 0. angleCosC(i,j,bi,bj) = 1. angleSinC(i,j,bi,bj) = 0. u2zonDir(i,j,bi,bj) = 1. v2zonDir(i,j,bi,bj) = 0. ENDDO cosFacU(j,bi,bj) = 1. cosFacV(j,bi,bj) = 1. sqCosFacU(j,bi,bj) = 1. sqCosFacV(j,bi,bj) = 1. ENDDO ENDDO ENDDO C Two examples are shown in this code. One illustrates the C initialization of a cartesian grid. The other shows the C inialization of a spherical polar grid. Other orthonormal grids C can be fitted into this design. In this case custom metric terms C also need adding to account for the projections of velocity C vectors onto these grids. The structure used here also makes it C possible to implement less regular grid mappings. In particular: C o Schemes which leave out blocks of the domain that are C all land could be supported. C o Multi-level schemes such as icosohedral or cubic C grid projectedions onto a sphere can also be fitted C within the strategy we use. C Both of the above also require modifying the support C routines that map computational blocks to simulation C domain blocks. C-- Set up horizontal grid and coordinate system IF ( usingCartesianGrid ) THEN CALL INI_CARTESIAN_GRID( myThid ) ELSEIF ( usingSphericalPolarGrid ) THEN CALL INI_SPHERICAL_POLAR_GRID( myThid ) ELSEIF ( usingCurvilinearGrid ) THEN CALL INI_CURVILINEAR_GRID( myThid ) ELSEIF ( usingCylindricalGrid ) THEN CALL INI_CYLINDER_GRID( myThid ) ELSE _BEGIN_MASTER(myThid) WRITE(msgBuf,'(2A)') 'S/R INI_GRID: ', & 'No grid coordinate system has been selected' CALL PRINT_ERROR( msgBuf , myThid) CALL ALL_PROC_DIE( 0 ) STOP 'ABNORMAL END: S/R INI_GRID' _END_MASTER(myThid) ENDIF C-- Calculate reciprocals grid lengths (formerly part of INI_MASKS_ETC) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF ( dxG(i,j,bi,bj) .NE. 0. ) & recip_dxG(i,j,bi,bj) = 1. _d 0/dxG(i,j,bi,bj) IF ( dyG(i,j,bi,bj) .NE. 0. ) & recip_dyG(i,j,bi,bj) = 1. _d 0/dyG(i,j,bi,bj) IF ( dxC(i,j,bi,bj) .NE. 0. ) & recip_dxC(i,j,bi,bj) = 1. _d 0/dxC(i,j,bi,bj) IF ( dyC(i,j,bi,bj) .NE. 0. ) & recip_dyC(i,j,bi,bj) = 1. _d 0/dyC(i,j,bi,bj) IF ( dxF(i,j,bi,bj) .NE. 0. ) & recip_dxF(i,j,bi,bj) = 1. _d 0/dxF(i,j,bi,bj) IF ( dyF(i,j,bi,bj) .NE. 0. ) & recip_dyF(i,j,bi,bj) = 1. _d 0/dyF(i,j,bi,bj) IF ( dxV(i,j,bi,bj) .NE. 0. ) & recip_dxV(i,j,bi,bj) = 1. _d 0/dxV(i,j,bi,bj) IF ( dyU(i,j,bi,bj) .NE. 0. ) & recip_dyU(i,j,bi,bj) = 1. _d 0/dyU(i,j,bi,bj) IF ( rA (i,j,bi,bj) .NE. 0. ) & recip_rA (i,j,bi,bj) = 1. _d 0/rA (i,j,bi,bj) IF ( rAs(i,j,bi,bj) .NE. 0. ) & recip_rAs(i,j,bi,bj) = 1. _d 0/rAs(i,j,bi,bj) IF ( rAw(i,j,bi,bj) .NE. 0. ) & recip_rAw(i,j,bi,bj) = 1. _d 0/rAw(i,j,bi,bj) IF ( rAz(i,j,bi,bj) .NE. 0. ) & recip_rAz(i,j,bi,bj) = 1. _d 0/rAz(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_MONITOR IF ( MASTER_CPU_IO(myThid) ) THEN C-- only the master thread is allowed to switch On/Off mon_write_stdout C & mon_write_mnc (since it is the only thread that uses those flags): IF (monitor_stdio) THEN mon_write_stdout = .TRUE. ELSE mon_write_stdout = .FALSE. ENDIF mon_write_mnc = .FALSE. #ifdef ALLOW_MNC IF (useMNC .AND. monitor_mnc) THEN DO i = 1,MAX_LEN_MBUF mon_fname(i:i) = ' ' ENDDO mon_fname(1:12) = 'monitor_grid' CALL MNC_CW_SET_UDIM(mon_fname, 1, myThid) mon_write_mnc = .TRUE. ENDIF #endif /* ALLOW_MNC */ ENDIF C Print out statistics of each horizontal grid array (helps when debugging) CALL MON_PRINTSTATS_RS(1,xC,'XC',myThid) CALL MON_PRINTSTATS_RS(1,xG,'XG',myThid) CALL MON_PRINTSTATS_RS(1,dxC,'DXC',myThid) CALL MON_PRINTSTATS_RS(1,dxF,'DXF',myThid) CALL MON_PRINTSTATS_RS(1,dxG,'DXG',myThid) CALL MON_PRINTSTATS_RS(1,dxV,'DXV',myThid) CALL MON_PRINTSTATS_RS(1,yC,'YC',myThid) CALL MON_PRINTSTATS_RS(1,yG,'YG',myThid) CALL MON_PRINTSTATS_RS(1,dyC,'DYC',myThid) CALL MON_PRINTSTATS_RS(1,dyF,'DYF',myThid) CALL MON_PRINTSTATS_RS(1,dyG,'DYG',myThid) CALL MON_PRINTSTATS_RS(1,dyU,'DYU',myThid) CALL MON_PRINTSTATS_RS(1,rA,'RA',myThid) CALL MON_PRINTSTATS_RS(1,rAw,'RAW',myThid) CALL MON_PRINTSTATS_RS(1,rAs,'RAS',myThid) CALL MON_PRINTSTATS_RS(1,rAz,'RAZ',myThid) CALL MON_PRINTSTATS_RS(1,angleCosC,'AngleCS',myThid) CALL MON_PRINTSTATS_RS(1,angleSinC,'AngleSN',myThid) IF ( MASTER_CPU_IO(myThid) ) THEN mon_write_stdout = .FALSE. mon_write_mnc = .FALSE. ENDIF #endif /* ALLOW_MONITOR */ C-- Everyone else must wait for the grid to be set _BARRIER RETURN END