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 CONSTANT
, only : OMEGA,RADIUS,TWOPI,SDAY,radian
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
DLAT=DLAT_DG*radian
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))
DXP(J) = RADIUS*DLON*COSP(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)
DXP1 = RADIUS*DLON*COSP1
DO J=J_0STG,J_1STG
COSV(J) = .5*(COSP(J-1)+COSP(J))
DXV(J) = .5*(DXP(J-1)+DXP(J))
DYV(J) = RADIUS*(LAT(J)-LAT(J-1))
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
DYP(1) = RADIUS*(LAT(2)-LAT(1)-0.5*DLAT)
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
DYP(JM) = RADIUS*(LAT(JM)-LAT(JM-1)-0.5*DLAT)
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**** adjacent velocity points
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