model/src/ini_grid.F

11033c8df1 Jean*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
924557e60a Chri*0003
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: INI_GRID
b6356366ca Ed H*0006
9366854e02 Chri*0007 C     !INTERFACE:
924557e60a Chri*0008       SUBROUTINE INI_GRID( myThid )
b6356366ca Ed H*0009 C     !DESCRIPTION:
                0010 C     These arrays are used throughout the code in evaluating gradients,
                0011 C     integrals and spatial avarages. This routine is called separately
                0012 C     by each thread and initializes only the region of the domain it is
                0013 C     "responsible" for.
924557e60a Chri*0014
d87b427b2d Jean*0015 C     !CALLING SEQUENCE:
                0016 C     INI_GRID
                0017 C      |   -- LOAD_GRID_SPACING
                0018 C      |   -- INI_VERTICAL_GRID
                0019 C      |    / INI_CARTESIAN_GRID
                0020 C      |   /  INI_SPHERICAL_POLAR_GRID
                0021 C      |   \  INI_CURVILINEAR_GRID
                0022 C      |    \ INI_CYLINDER_GRID
                0023
9366854e02 Chri*0024 C     !USES:
                0025       IMPLICIT NONE
924557e60a Chri*0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
a9e63588d9 Jean*0029 #include "GRID.h"
f31930e56f Ed H*0030 #ifdef ALLOW_MNC
                0031 #include "MNC_PARAMS.h"
                0032 #endif
b6356366ca Ed H*0033 #ifdef ALLOW_MONITOR
                0034 #include "MONITOR.h"
                0035 #endif
924557e60a Chri*0036
9366854e02 Chri*0037 C     !INPUT/OUTPUT PARAMETERS:
b9306711eb Jean*0038 C     myThid  :: my Thread Id number
924557e60a Chri*0039       INTEGER myThid
                0040
67eed25d83 Jean*0041 C     !FUNCTIONS:
49aab2cab9 Jean*0042       LOGICAL  MASTER_CPU_IO
                0043       EXTERNAL MASTER_CPU_IO
                0044
9366854e02 Chri*0045 C     !LOCAL VARIABLES:
67eed25d83 Jean*0046 C     bi, bj  :: tile indices
                0047 C     i, j    :: Loop counters
b9306711eb Jean*0048 C     msgBuf  :: Informational/error message buffer
67eed25d83 Jean*0049       INTEGER bi, bj
                0050       INTEGER i, j
924557e60a Chri*0051       CHARACTER*(MAX_LEN_MBUF) msgBuf
67eed25d83 Jean*0052 CEOP
                0053
                0054 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
924557e60a Chri*0055
71b5d4795f Jean*0056 C--   Load grid spacing (vector) from files
f76021331f Jean*0057       CALL LOAD_GRID_SPACING( myThid )
                0058
71b5d4795f Jean*0059 C--   Set up vertical grid and coordinate system
a74b31ee34 Alis*0060       CALL INI_VERTICAL_GRID( myThid )
                0061
71b5d4795f Jean*0062 C--   Initialise (everywhere) all horizontal grid array to null value
                0063 C     Note: some arrays are not defined in some parts of the halo
                0064 C     region. We set them to zero here for safety. If they are ever
                0065 C     referred to, especially in the denominator then it is a mistake!
                0066       DO bj = myByLo(myThid), myByHi(myThid)
                0067        DO bi = myBxLo(myThid), myBxHi(myThid)
                0068         DO j=1-OLy,sNy+OLy
                0069          DO i=1-OLx,sNx+OLx
                0070           xC(i,j,bi,bj)  = 0.
                0071           yC(i,j,bi,bj)  = 0.
                0072           xG(i,j,bi,bj)  = 0.
                0073           yG(i,j,bi,bj)  = 0.
                0074           dxC(i,j,bi,bj) = 0.
                0075           dyC(i,j,bi,bj) = 0.
                0076           dxG(i,j,bi,bj) = 0.
                0077           dyG(i,j,bi,bj) = 0.
                0078           dxF(i,j,bi,bj) = 0.
                0079           dyF(i,j,bi,bj) = 0.
                0080           dxV(i,j,bi,bj) = 0.
                0081           dyU(i,j,bi,bj) = 0.
                0082           rA(i,j,bi,bj)  = 0.
                0083           rAz(i,j,bi,bj) = 0.
                0084           rAw(i,j,bi,bj) = 0.
                0085           rAs(i,j,bi,bj) = 0.
                0086           recip_dxG(i,j,bi,bj) = 0.
                0087           recip_dyG(i,j,bi,bj) = 0.
                0088           recip_dxC(i,j,bi,bj) = 0.
                0089           recip_dyC(i,j,bi,bj) = 0.
                0090           recip_dxF(i,j,bi,bj) = 0.
                0091           recip_dyF(i,j,bi,bj) = 0.
                0092           recip_dxV(i,j,bi,bj) = 0.
                0093           recip_dyU(i,j,bi,bj) = 0.
                0094           recip_rA (i,j,bi,bj) = 0.
                0095           recip_rAs(i,j,bi,bj) = 0.
                0096           recip_rAw(i,j,bi,bj) = 0.
                0097           recip_rAz(i,j,bi,bj) = 0.
                0098           tanPhiAtU(i,j,bi,bj) = 0.
                0099           tanPhiAtV(i,j,bi,bj) = 0.
                0100           angleCosC(i,j,bi,bj) = 1.
                0101           angleSinC(i,j,bi,bj) = 0.
                0102           u2zonDir(i,j,bi,bj)  = 1.
                0103           v2zonDir(i,j,bi,bj)  = 0.
                0104          ENDDO
                0105          cosFacU(j,bi,bj)   = 1.
                0106          cosFacV(j,bi,bj)   = 1.
                0107          sqCosFacU(j,bi,bj) = 1.
                0108          sqCosFacV(j,bi,bj) = 1.
                0109         ENDDO
                0110        ENDDO
                0111       ENDDO
                0112
b6356366ca Ed H*0113 C     Two examples are shown in this code. One illustrates the
                0114 C     initialization of a cartesian grid. The other shows the
                0115 C     inialization of a spherical polar grid. Other orthonormal grids
                0116 C     can be fitted into this design. In this case custom metric terms
                0117 C     also need adding to account for the projections of velocity
                0118 C     vectors onto these grids.  The structure used here also makes it
                0119 C     possible to implement less regular grid mappings. In particular:
b9306711eb Jean*0120 C      o Schemes which leave out blocks of the domain that are
                0121 C        all land could be supported.
                0122 C      o Multi-level schemes such as icosohedral or cubic
                0123 C        grid projectedions onto a sphere can also be fitted
                0124 C       within the strategy we use.
                0125 C        Both of the above also require modifying the support
                0126 C        routines that map computational blocks to simulation
                0127 C        domain blocks.
b6356366ca Ed H*0128
71b5d4795f Jean*0129 C--   Set up horizontal grid and coordinate system
924557e60a Chri*0130       IF ( usingCartesianGrid ) THEN
b6356366ca Ed H*0131         CALL INI_CARTESIAN_GRID( myThid )
924557e60a Chri*0132       ELSEIF ( usingSphericalPolarGrid ) THEN
b6356366ca Ed H*0133         CALL INI_SPHERICAL_POLAR_GRID( myThid )
aea29c8517 Alis*0134       ELSEIF ( usingCurvilinearGrid ) THEN
b6356366ca Ed H*0135         CALL INI_CURVILINEAR_GRID( myThid )
0ac260a803 Andr*0136       ELSEIF ( usingCylindricalGrid ) THEN
b6aef9f254 Jean*0137         CALL INI_CYLINDER_GRID( myThid )
924557e60a Chri*0138       ELSE
b6356366ca Ed H*0139         _BEGIN_MASTER(myThid)
                0140         WRITE(msgBuf,'(2A)') 'S/R INI_GRID: ',
                0141      &       'No grid coordinate system has been selected'
                0142         CALL PRINT_ERROR( msgBuf , myThid)
67eed25d83 Jean*0143         CALL ALL_PROC_DIE( 0 )
b6356366ca Ed H*0144         STOP 'ABNORMAL END: S/R INI_GRID'
                0145         _END_MASTER(myThid)
924557e60a Chri*0146       ENDIF
aea29c8517 Alis*0147
71b5d4795f Jean*0148 C--   Calculate reciprocals grid lengths (formerly part of INI_MASKS_ETC)
67eed25d83 Jean*0149       DO bj = myByLo(myThid), myByHi(myThid)
                0150        DO bi = myBxLo(myThid), myBxHi(myThid)
                0151         DO j=1-OLy,sNy+OLy
                0152          DO i=1-OLx,sNx+OLx
                0153           IF ( dxG(i,j,bi,bj) .NE. 0. )
                0154      &    recip_dxG(i,j,bi,bj) = 1. _d 0/dxG(i,j,bi,bj)
                0155           IF ( dyG(i,j,bi,bj) .NE. 0. )
                0156      &    recip_dyG(i,j,bi,bj) = 1. _d 0/dyG(i,j,bi,bj)
                0157           IF ( dxC(i,j,bi,bj) .NE. 0. )
                0158      &    recip_dxC(i,j,bi,bj) = 1. _d 0/dxC(i,j,bi,bj)
                0159           IF ( dyC(i,j,bi,bj) .NE. 0. )
                0160      &    recip_dyC(i,j,bi,bj) = 1. _d 0/dyC(i,j,bi,bj)
                0161           IF ( dxF(i,j,bi,bj) .NE. 0. )
                0162      &    recip_dxF(i,j,bi,bj) = 1. _d 0/dxF(i,j,bi,bj)
                0163           IF ( dyF(i,j,bi,bj) .NE. 0. )
                0164      &    recip_dyF(i,j,bi,bj) = 1. _d 0/dyF(i,j,bi,bj)
                0165           IF ( dxV(i,j,bi,bj) .NE. 0. )
                0166      &    recip_dxV(i,j,bi,bj) = 1. _d 0/dxV(i,j,bi,bj)
                0167           IF ( dyU(i,j,bi,bj) .NE. 0. )
                0168      &    recip_dyU(i,j,bi,bj) = 1. _d 0/dyU(i,j,bi,bj)
                0169           IF ( rA (i,j,bi,bj) .NE. 0. )
                0170      &    recip_rA (i,j,bi,bj) = 1. _d 0/rA (i,j,bi,bj)
                0171           IF ( rAs(i,j,bi,bj) .NE. 0. )
                0172      &    recip_rAs(i,j,bi,bj) = 1. _d 0/rAs(i,j,bi,bj)
                0173           IF ( rAw(i,j,bi,bj) .NE. 0. )
                0174      &    recip_rAw(i,j,bi,bj) = 1. _d 0/rAw(i,j,bi,bj)
                0175           IF ( rAz(i,j,bi,bj) .NE. 0. )
                0176      &    recip_rAz(i,j,bi,bj) = 1. _d 0/rAz(i,j,bi,bj)
                0177          ENDDO
                0178         ENDDO
                0179        ENDDO
                0180       ENDDO
                0181
                0182 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0183
a3159bd164 Patr*0184 #ifdef ALLOW_MONITOR
49aab2cab9 Jean*0185       IF ( MASTER_CPU_IO(myThid) ) THEN
a9e63588d9 Jean*0186 C--   only the master thread is allowed to switch On/Off mon_write_stdout
ff02675122 Jean*0187 C     & mon_write_mnc (since it is the only thread that uses those flags):
a9e63588d9 Jean*0188
                0189         IF (monitor_stdio) THEN
                0190           mon_write_stdout = .TRUE.
                0191         ELSE
                0192           mon_write_stdout = .FALSE.
                0193         ENDIF
                0194         mon_write_mnc = .FALSE.
b6356366ca Ed H*0195 #ifdef ALLOW_MNC
a9e63588d9 Jean*0196         IF (useMNC .AND. monitor_mnc) THEN
                0197           DO i = 1,MAX_LEN_MBUF
                0198             mon_fname(i:i) = ' '
                0199           ENDDO
                0200           mon_fname(1:12) = 'monitor_grid'
                0201           CALL MNC_CW_SET_UDIM(mon_fname, 1, myThid)
                0202           mon_write_mnc = .TRUE.
                0203         ENDIF
b6356366ca Ed H*0204 #endif /*  ALLOW_MNC  */
a9e63588d9 Jean*0205
                0206       ENDIF
                0207
67eed25d83 Jean*0208 C     Print out statistics of each horizontal grid array (helps when debugging)
a9e63588d9 Jean*0209       CALL MON_PRINTSTATS_RS(1,xC,'XC',myThid)
                0210       CALL MON_PRINTSTATS_RS(1,xG,'XG',myThid)
                0211       CALL MON_PRINTSTATS_RS(1,dxC,'DXC',myThid)
                0212       CALL MON_PRINTSTATS_RS(1,dxF,'DXF',myThid)
                0213       CALL MON_PRINTSTATS_RS(1,dxG,'DXG',myThid)
                0214       CALL MON_PRINTSTATS_RS(1,dxV,'DXV',myThid)
                0215       CALL MON_PRINTSTATS_RS(1,yC,'YC',myThid)
                0216       CALL MON_PRINTSTATS_RS(1,yG,'YG',myThid)
                0217       CALL MON_PRINTSTATS_RS(1,dyC,'DYC',myThid)
                0218       CALL MON_PRINTSTATS_RS(1,dyF,'DYF',myThid)
                0219       CALL MON_PRINTSTATS_RS(1,dyG,'DYG',myThid)
                0220       CALL MON_PRINTSTATS_RS(1,dyU,'DYU',myThid)
                0221       CALL MON_PRINTSTATS_RS(1,rA,'RA',myThid)
                0222       CALL MON_PRINTSTATS_RS(1,rAw,'RAW',myThid)
                0223       CALL MON_PRINTSTATS_RS(1,rAs,'RAS',myThid)
                0224       CALL MON_PRINTSTATS_RS(1,rAz,'RAZ',myThid)
b6aef9f254 Jean*0225       CALL MON_PRINTSTATS_RS(1,angleCosC,'AngleCS',myThid)
                0226       CALL MON_PRINTSTATS_RS(1,angleSinC,'AngleSN',myThid)
b6356366ca Ed H*0227
49aab2cab9 Jean*0228       IF ( MASTER_CPU_IO(myThid) ) THEN
a9e63588d9 Jean*0229         mon_write_stdout = .FALSE.
                0230         mon_write_mnc    = .FALSE.
                0231       ENDIF
                0232 #endif /* ALLOW_MONITOR */
49aab2cab9 Jean*0233
                0234 C--   Everyone else must wait for the grid to be set
                0235       _BARRIER
398b491097 Alis*0236
924557e60a Chri*0237       RETURN
                0238       END