```

MODULE GEOM 106,2
!@sum  GEOM contains spherical geometric variables and arrays
!@auth Original development team
!@ver  1.0 (B grid version)
!@cont GEOM_B
USE MODEL_COM, only : IM,JM,LM,FIM,BYIM
!     USE DOMAIN_DECOMP, only : grid, get
!     USE DOMAIN_DECOMP, only : halo_update, north, south, checksum
IMPLICIT NONE
C**** The primary grid is the A grid (including both poles)
C**** The secondary grid is for the B grid velocities, located on the
C**** vertices of the A grid (Note: first velocity point is J=2)
C**** Polar boxes can have different latitudinal size and are treated
C**** as though they were 1/IM of their actual area
SAVE
!@param  DLON grid spacing in longitude (deg)
REAL*8, PARAMETER :: DLON = TWOPI*BYIM
C**** For the wonderland model set DLON=DLON/3
c      REAL*8, PARAMETER :: DLON=TWOPI/(IM*3)
!@param  DLAT,DLAT_DG grid spacing in latitude (rad,deg)
REAL*8  :: DLAT,DLAT_DG
!@param  FJEQ equatorial value of J
REAL*8, PARAMETER :: FJEQ=.5*(1+JM)
!@var  J1U index of southernmost latitude (currently 2, later 1)
INTEGER, parameter :: J1U = 2
!@var  JRANGE_HEMI lowest,highest lat index for SH,NH for A,B grid
INTEGER, parameter, dimension(2,2,2) :: JRANGE_HEMI = reshape(
*  (/1,JM/2,  1+JM/2,JM,  J1U,J1U-1+JM/2, J1U-1+JM/2,JM+J1U-2/),
*  (/2,2,2/))
!@var  LAT latitude of mid point of primary grid box (radians)
REAL*8, ALLOCATABLE, DIMENSION(:) :: LAT
!@var  LAT_DG latitude of mid points of primary and sec. grid boxs (deg)
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: LAT_DG
!@var  LON longitude of mid points of primary grid box (radians)
REAL*8, ALLOCATABLE, DIMENSION(:) :: LON
!@var  LON_DG longitude of mid points of prim. and sec. grid boxes (deg)
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: LON_DG
!@var  DXYP,BYDXYP area of grid box (+inverse) (m^2)
C**** Note that this is not the exact area, but is what is required for
C**** some B-grid conservation quantities
C**** Decomposition exception: DXYP not distributed because individual elements
C**** of the 1D array are needed by all processes.
REAL*8  DXYP(JM)
REAL*8, ALLOCATABLE, DIMENSION(:) :: BYDXYP
!@var AREAG global integral of area (m^2)
REAL*8 :: AREAG
!@var WTJ area weighting used in JLMAP, JKMAP (for hemispheric means)
REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: WTJ

!@var DXYV,BYDXYV area of grid box around velocity point (recip.)(m^2)
REAL*8, ALLOCATABLE, DIMENSION(:) ::
&                  DXYV,BYDXYV

!@var  DXP,DYP,BYDXP,BYDYP distance between points on primary grid
!@+     (+inverse)
REAL*8, ALLOCATABLE, DIMENSION(:) ::
&                  DXP,DYP,BYDXP,BYDYP
!@var  DXV,DYV distance between velocity points (secondary grid)
REAL*8, ALLOCATABLE, DIMENSION(:) :: DXV,DYV
!@var  DXYN,DXYS half box areas to the North,South of primary grid point
REAL*8, ALLOCATABLE, DIMENSION(:) :: DXYS,DXYN
!@var  SINP sin of latitude at primary grid points
REAL*8, ALLOCATABLE, DIMENSION(:) :: SINP
!@var  COSP, COSV cos of latitude at primary, secondary latitudes
REAL*8, ALLOCATABLE, DIMENSION(:) :: COSP,COSV
!@var  RAPVS,RAPVN,RAVPS,RAVPN area scalings for primary and sec. grid
REAL*8, ALLOCATABLE, DIMENSION(:) ::
&                  RAPVS,RAPVN,RAVPS,RAVPN

!@var SINIV,COSIV,SINIP,COSIP longitud. sin,cos for wind,pressure grid
REAL*8, ALLOCATABLE, DIMENSION(:) :: SINIV,COSIV,SINIP,COSIP
!@var  RAVJ scaling for A grid U/V to B grid points (func. of lat. j)
!@var  RAPJ scaling for B grid -> A grid conversion (1/4,1/im at poles)
REAL*8, ALLOCATABLE, DIMENSION(:,:) ::
&                  RAPJ,RAVJ
!@var  IDJJ J index of adjacent U/V points for A grid (func. of lat. j)
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: IDJJ
!@var  IDIJ I index of adjacent U/V points for A grid (func. of lat/lon)
INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::
&                  IDIJ
!@var  KMAXJ varying number of adjacent velocity points
INTEGER, ALLOCATABLE, DIMENSION(:) :: KMAXJ
!@var  IMAXJ varying number of used longitudes
INTEGER, ALLOCATABLE, DIMENSION(:) :: IMAXJ
!@var  FCOR latitudinally varying coriolis parameter
REAL*8, ALLOCATABLE, DIMENSION(:) :: FCOR

CONTAINS

SUBROUTINE GEOM_B 1,3
!@sum  GEOM_B Calculate spherical geometry for B grid
!@auth Original development team (modifications by G. Schmidt)
!@ver  1.0 (B grid version)
USE DOMAIN_DECOMP, only : grid, get
USE DOMAIN_DECOMP, only : halo_update, north, south, checksum
IMPLICIT NONE
REAL*8, PARAMETER :: EDPERD=1.,EDPERY = 365.

INTEGER :: I,J,K,IM1  !@var I,J,K,IM1  loop variables
INTEGER :: JVPO,JMHALF
REAL*8  :: RAVPO,LAT1,COSP1,DXP1
INTEGER J_0,J_1,J_0S,J_1S,J_0STG,J_1STG

C**** latitudinal spacing depends on whether you have even spacing or
C**** a partial box at the pole
DLAT_DG=180./JM                     ! even spacing (default)
IF (JM.eq.46) DLAT_DG=180./(JM-1)   ! 1/2 box at pole for 4x5
cc    IF (JM.eq.24) DLAT_DG=180./(JM-1)   ! 1/2 box at pole, orig 8x10
IF (JM.eq.24) DLAT_DG=180./(JM-1.5) ! 1/4 box at pole, 'real' 8x10
CALL GET(grid, J_STRT     =J_0,    J_STOP     =J_1,
&               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S,
&               J_STRT_STGR=J_0STG, j_STOP_STGR=J_1STG)

IF (grid%HAVE_SOUTH_POLE) THEN
LAT(1)  = -.25*TWOPI
SINP(1)  = -1.
COSP(1)  = 0.
DXP(1)  = 0.
END IF
IF (grid%HAVE_NORTH_POLE) THEN
LAT(JM) = .25*TWOPI
SINP(JM) = 1.
COSP(JM) = 0.
DXP(JM) = 0.
END IF
DO J=J_0S,J_1S
LAT(J)  = DLAT*(J-FJEQ)
SINP(J) = SIN(LAT(J))
COSP(J) = COS(LAT(J))
END DO
BYDXP(J_0S:J_1S) = 1.D0/DXP(J_0S:J_1S)
C****Update halos for arrays cosp, dxp, and lat
CALL CHECKSUM(grid, COSP, 136, "GEOM_B.f")
CALL CHECKSUM(grid, DXP , 137, "GEOM_B.f")
CALL CHECKSUM(grid, LAT , 138, "GEOM_B.f")
CALL HALO_UPDATE( grid, COSP, from=SOUTH)
CALL HALO_UPDATE( grid, DXP , from=SOUTH)
CALL HALO_UPDATE( grid, LAT , from=SOUTH)

LAT1    = DLAT*(1.-FJEQ)
COSP1   = COS(LAT1)

DO J=J_0STG,J_1STG
COSV(J) = .5*(COSP(J-1)+COSP(J))
DXV(J)  = .5*(DXP(J-1)+DXP(J))
C**** The following corrections have no effect for half polar boxes
C**** but are important for full and quarter polar box cases.
IF (J.eq.2) THEN
COSV(J) = .5*(COSP1+COSP(J))
DXV(J)  = .5*(DXP1+DXP(J))
END IF
IF (J.eq.JM) THEN
COSV(J) = .5*(COSP(J-1)+COSP1)
DXV(J)  = .5*(DXP(J-1)+DXP1)
END IF
C****
END DO
C****POLES
IF (grid%HAVE_SOUTH_POLE) THEN
DXYP(1) = .5*DXV(2)*DYP(1)
BYDXYP(1) = 1./DXYP(1)
DXYS(1)  = 0.
DXYN(1)  = DXYP(1)
ENDIF
IF (grid%HAVE_NORTH_POLE) THEN
DXYP(JM)= .5*DXV(JM)*DYP(JM)
BYDXYP(JM) = 1./DXYP(JM)
DXYS(JM) = DXYP(JM)
DXYN(JM) = 0.
END IF
AREAG = DXYP(1)+DXYP(JM)

DO J=J_0S,J_1S
DYP(J)  = .5*(DYV(J)+DYV(J+1))
DXYP(J) = .5*(DXV(J)+DXV(J+1))*DYP(J)
BYDXYP(J) = 1./DXYP(J)
DXYS(J) = .5*DXYP(J)
DXYN(J) = .5*DXYP(J)
AREAG = AREAG+DXYP(J)
END DO
BYDYP(:) = 1.D0/DYP(:)
AREAG = AREAG*FIM
C****EXCEPTION!
C****      dxyp is not distributed -> make sure all procs. have right values for
C****                                 the whole array.
C****  CALL ALL_GATHER(DXYP)

C****POLES
IF (grid%HAVE_SOUTH_POLE) THEN
RAVPS(1)  = 0.
RAPVS(1)  = 0.
END IF
IF (grid%HAVE_NORTH_POLE) THEN
RAVPN(JM) = 0.
RAPVN(JM) = 0.
END IF
DO J=J_0STG,J_1STG
DXYV(J) = DXYN(J-1)+DXYS(J)
BYDXYV(J) = 1./DXYV(J)
RAPVS(J)   = .5*DXYS(J)/DXYV(J)
RAVPS(J)   = .5*DXYS(J)/DXYP(J)
END DO
CALL CHECKSUM(grid,RAPVN, 210, "GEOM_B.f")
CALL HALO_UPDATE(grid,RAPVN, from=NORTH)
DO J=J_0,J_1S
RAPVN(J) = .5*DXYN(J)/DXYV(J+1)
RAVPN(J) = .5*DXYN(J)/DXYP(J)
END DO
C**** LONGITUDES (degrees); used in ILMAP
LON_DG(1,1) = -180.+360./(2.*FLOAT(IM))
LON_DG(1,2) = -180.+360./    FLOAT(IM)
DO I=2,IM
LON_DG(I,1) = LON_DG(I-1,1)+360./FLOAT(IM)
LON_DG(I,2) = LON_DG(I-1,2)+360./FLOAT(IM)
END DO
C**** LATITUDES (degrees); used extensively in the diagn. print routines
IF (grid%HAVE_SOUTH_POLE) LAT_DG(1,1:2)=-90.
IF (grid%HAVE_NORTH_POLE) LAT_DG(JM,1)=90.
DO J=J_0S,J_1S
LAT_DG(J,1)=DLAT_DG*(J-FJEQ)    ! primary (tracer) latitudes
END DO
DO J=J_0STG,J_1STG
LAT_DG(J,2)=DLAT_DG*(J-JM/2-1)  ! secondary (velocity) latitudes
END DO
C**** WTJ: area weighting for JKMAP, JLMAP hemispheres
JMHALF= JM/2
DO J=J_0,J_1
WTJ(J,1,1)=1.
WTJ(J,2,1)=2.*FIM*DXYP(J)/AREAG
END DO
DO J=J_0STG,J_1STG
WTJ(J,1,2)=1.
WTJ(J,2,2)=2.*FIM*DXYV(J)/AREAG
END DO
C****EQUATOR
IF (grid%HAVE_EQUATOR) THEN
WTJ(JMHALF+1,1,2)=.5
WTJ(JMHALF+1,2,2)=WTJ(JMHALF+1,2,2)/2.
END IF
C****POLE
IF (grid%HAVE_SOUTH_POLE) THEN
WTJ(1,1,2)=0.
WTJ(1,2,2)=0.
END IF

C**** CALCULATE CORIOLIS PARAMETER (NOW RESOLUTION INDEPENDENT)
c      OMEGA = TWOPI*(EDPERD+EDPERY)/(EDPERD*EDPERY*SDAY)
IF (grid%HAVE_SOUTH_POLE)
&      FCOR(1)  = -OMEGA*DXV(2)*DXV(2)/DLON
IF (grid%HAVE_NORTH_POLE)
&      FCOR(JM) = OMEGA*DXV(JM)*DXV(JM)/DLON
C****Update halo for DXV array
CALL CHECKSUM(grid,DXV, 260, "GEOM_B.f")
CALL HALO_UPDATE(grid,DXV, from=NORTH)
DO J=J_0S,J_1S
FCOR(J) = OMEGA*(DXV(J)*DXV(J)-DXV(J+1)*DXV(J+1))/DLON
END DO

C**** Set indexes and scalings for the influence of A grid points on

C**** Calculate relative directions of polar box to nearby U,V points
DO I=1,IM
SINIV(I)=SIN((I-1)*DLON)
COSIV(I)=COS((I-1)*TWOPI*BYIM) ! DLON)
LON(I)=DLON*(I-.5)
SINIP(I)=SIN(LON(I))
COSIP(I)=COS(LON(I))
END DO

C**** Conditions at the poles
DO J=1,JM,JM-1
IF (((J .EQ. 1) .AND. (grid%HAVE_SOUTH_POLE)) .OR.
&      ((J .EQ.JM) .AND. (grid%HAVE_NORTH_POLE))) THEN
IF(J.EQ.1) THEN
JVPO=2
RAVPO=2.*RAPVN(1)
ELSE
JVPO=JM
RAVPO=2.*RAPVS(JM)
END IF
KMAXJ(J)=IM
IMAXJ(J)=1
RAVJ(1:KMAXJ(J),J)=RAVPO
RAPJ(1:KMAXJ(J),J)=BYIM
IDJJ(1:KMAXJ(J),J)=JVPO
DO K=1,KMAXJ(J)
IDIJ(K,1:IM,J)=K
END DO
END IF
END DO
C**** Conditions at non-polar points
DO J=J_0S,J_1S
KMAXJ(J)=4
IMAXJ(J)=IM
DO K=1,2
RAVJ(K,J)=RAPVS(J)
RAPJ(K,J)=RAVPS(J)    ! = .25
IDJJ(K,J)=J
RAVJ(K+2,J)=RAPVN(J)
RAPJ(K+2,J)=RAVPN(J)  ! = .25
IDJJ(K+2,J)=J+1
END DO
IM1=IM
DO I=1,IM
IDIJ(1,I,J)=IM1
IDIJ(2,I,J)=I
IDIJ(3,I,J)=IM1
IDIJ(4,I,J)=I
IM1=I
END DO
END DO

RETURN
END SUBROUTINE GEOM_B

END MODULE GEOM

SUBROUTINE ALLOC_GEOM(grd_dum) 1,4
USE MODEL_COM, only : IM
USE DOMAIN_DECOMP, only: DYN_GRID, get
USE GEOM, only: LAT,LAT_DG,
&                LON,LON_DG,
&                DXYP,BYDXYP,WTJ,
&                DXYV,BYDXYV,
&                DXP,DYP,BYDXP,BYDYP,
&                DXV,DYV,DXYS,DXYN,
&                SINP,COSP,COSV,
&                RAPVS,RAPVN,RAVPS,RAVPN,
&                SINIV,COSIV,SINIP,COSIP,
&                RAPJ,RAVJ,
&                IDJJ,IDIJ,KMAXJ,IMAXJ,
&                FCOR

TYPE (DYN_GRID),   INTENT(IN)    :: grd_dum

CALL GET(grd_dum,J_STRT_HALO=J_0H,
&                 J_STOP_HALO=J_1H)

ALLOCATE( LAT       (J_0H:J_1H),
&          LAT_DG    (J_0H:J_1H, 2),
&          LON    (IM),
&          LON_DG (IM,           2),
&                                    STAT = IER)
ALLOCATE(
&          BYDXYP    (J_0H:J_1H),
&          WTJ       (J_0H:J_1H, 2, 2),
&          DXYV      (J_0H:J_1H),
&          BYDXYV    (J_0H:J_1H),
&          DXP       (J_0H:J_1H),
&          DYP       (J_0H:J_1H),
&          BYDXP     (J_0H:J_1H),
&          BYDYP     (J_0H:J_1H),
&          DXV       (J_0H:J_1H),
&          DYV       (J_0H:J_1H),
&          DXYS      (J_0H:J_1H),
&          DXYN      (J_0H:J_1H),
&          SINP      (J_0H:J_1H),
&          COSP      (J_0H:J_1H),
&          COSV      (J_0H:J_1H),
&          RAPVS     (J_0H:J_1H),
&          RAPVN     (J_0H:J_1H),
&          RAVPS     (J_0H:J_1H),
&          RAVPN     (J_0H:J_1H),
&                                    STAT = IER)
ALLOCATE(
&          SINIV  (IM),
&          COSIV  (IM),
&          SINIP  (IM),
&          COSIP  (IM),
&          RAPJ   (IM,J_0H:J_1H),
&          RAVJ   (IM,J_0H:J_1H),
&          IDJJ   (IM,J_0H:J_1H),
&          IDIJ(IM,IM,J_0H:J_1H),
&          KMAXJ     (J_0H:J_1H),
&          IMAXJ     (J_0H:J_1H),
&          FCOR      (J_0H:J_1H),
&                                    STAT = IER)

RETURN
END SUBROUTINE ALLOC_GEOM
```