suleg_mod.F90

Go to the documentation of this file.

MODULE suleg_mod
CONTAINS
SUBROUTINE suleg

USE parkind1  ,ONLY : jpim     ,jprb
USE parkind2  ,ONLY : jprh

USE tpm_gen
USE tpm_dim
USE tpm_constants
USE tpm_distr
USE tpm_fields

!USE SUGAW_MOD
USE supol_mod
USE sutrle_mod

#ifdef DOC

!**** *SULEG * - initialize the Legendre polynomials

!     Purpose.
!     --------
!           Initialize COMMON YOMLEG

!**   Interface.
!     ----------
!        *CALL* *SULEG*

!        Explicit arguments :
!        --------------------

!        Implicit arguments :
!        --------------------
!              COMMON YOMLEG

!     Method.
!     -------
!        See documentation

!     Externals.
!     ----------
!        SUGAW  (Gaussian latitudes)
!        SUPOLM (polynomials)
!        LFI routines for external IO's
!        Called by SUGEM.

!     Reference.
!     ----------
!        ECMWF Research Department documentation of the IFS

!     Author.
!     -------
!        Mats Hamrud and Philippe Courtier  *ECMWF*

!     Modifications.
!     --------------
!        Original : 87-10-15
!        MODIFICATION : 91-04  J.M. Piriou:
!                       - Read gaussian latitudes and PNM on LFI
!                       - If file missing, computes
!                       91-04 M.Hamrud:
!                       - IO Scheme introduced
!        MODIFICATION : 91-07-03 P.Courtier suppress derivatives
!        MODIFICATION : 91-07-03 P.Courtier computes RATATH and RACTHE
!        MODIFICATION : 91-07-03 P.Courtier change upper limit (NSMAX+1)
!        MODIFICATION : 91-07-03 P.Courtier change ordering
!     Order of the PNM in the file, as in the model :
!         - increasing wave numbers m
!         - for a given m, from n=NSMAX+1 to m
!        MODIFICATION : 92-07-02 R. Bubnova: shift RATATH calculation
!                                            to SUGEM1
!        MODIFICATION : 92-12-17 P.Courtier multitask computations
!        Modified by R. EL Khatib : 93-04-02 Set-up defaults controled by LECMWF
!        MODIFICATION : 93-03-19 D.Giard : n <= NTMAX
!        K. YESSAD    : 93-05-11 : DLMU --> global array DRMU(NDGSA:NDGEN).
!                       (not stored currently on LFI files).
!        MODIFICATION : 94-02-03 R. El Khatib : subroutine SULEG2 to write out
!                       the Leg. polynomials on workfile or LFI file
!        Modification : 94-08-31 M. Tolstykh: Setup for CUD interpolation
!        Modified by K. YESSAD (MARCH 1995): Extra-latitudes computations
!                 according to value of NDGSUR and LRPOLE only.
!                 + change fancy loop numbering.
!        Modified 98-08-10 by K. YESSAD: removal of LRPOLE option.
!          - removal of LRPOLE in YOMCT0.
!          - removal of code under LRPOLE.
!     ------------------------------------------------------------------

#endif


IMPLICIT NONE


!     ------------------------------------------------------------------
REAL(KIND=JPRB),ALLOCATABLE :: ZPNMG(:,:)

REAL(KIND=JPRH) :: DLRMU(r%ndgl)
REAL(KIND=JPRH) :: DLC(0:r%ntmax+1,0:r%ntmax+1)
REAL(KIND=JPRH) :: DLD(0:r%ntmax+1,0:r%ntmax+1)
REAL(KIND=JPRH) :: DLE(0:r%ntmax+1,0:r%ntmax+1)
REAL(KIND=JPRH) :: DLA(0:r%ntmax+1),DLB(0:r%ntmax+1),DLF(0:r%ntmax+1)
REAL(KIND=JPRH) :: DLG(0:r%ntmax+1),DLH(0:r%ntmax+1),DLI(0:r%ntmax+1)
REAL(KIND=JPRH) :: DLPOL(0:r%ntmax+1,0:r%ntmax+1)
!     ------------------------------------------------------------------

INTEGER(KIND=JPIM), PARAMETER :: JPKS=kind(zpnmg)
INTEGER(KIND=JPIM), PARAMETER :: JPKD=kind(dlg)

!     ------------------------------------------------------------------
REAL(KIND=JPRH) :: DA,DC,DD,DE
INTEGER(KIND=JPIM) :: KKN, KKM

!     LOCAL 
INTEGER(KIND=JPIM) :: IGLLOC, INM, IM , ICOUNT,&
             &JGL,  JM, JMLOC, JN, JNM


LOGICAL :: LLP1,LLP2


dc(kkn,kkm)=sqrt( (REAL(2*kkn+1,jpkd)*REAL(kkn+kkm-1,jpkd)&
                   &*REAL(KKN+KKM-3,JPKD))&
                &/ (REAL(2*KKN-3,JPKD)*REAL(KKN+KKM,JPKD)&
                   &*REAL(KKN+KKM-2,JPKD)) )
DD(kkn,kkm)=sqrt( (REAL(2*kkn+1,jpkd)*REAL(kkn+kkm-1,jpkd)&
                   &*REAL(KKN-KKM+1,JPKD))&
                &/ (REAL(2*KKN-1,JPKD)*REAL(KKN+KKM,JPKD)&
                  &*REAL(KKN+KKM-2,JPKD)) )
DE(kkn,kkm)=sqrt( (REAL(2*kkn+1,jpkd)*REAL(kkn-kkm,jpkd))&
                &/ (REAL(2*KKN-1,JPKD)*REAL(KKN+KKM,JPKD)) )
DA(kkn,kkm)=sqrt( (REAL(2*kkn+1,jpkd)*REAL(kkn-kkm,jpkd)&
                   &*REAL(KKN+KKM,JPKD))&
                &/  REAL(2*KKN-1,JPKD) )

!     ------------------------------------------------------------------
ALLOCATE(zpnmg(r%NSPOLEG,d%NLEI3D))

!*       0.    Some initializations.
!              ---------------------

llp1 = nprintlev>0
llp2 = nprintlev>1
IF(llp1) WRITE(nout,*) '=== ENTER ROUTINE SULEG ==='

!CALL GSTATS(140,0)  !MPL 4.12.08
ALLOCATE(f%RPNM(r%NLEI3,d%NSPOLEGL))
IF (llp2) WRITE(nout,9) 'F%RPNM    ',SIZE(f%RPNM),shape(f%RPNM) 
ALLOCATE(f%RMU(r%NDGL))
IF (llp2) WRITE(nout,9) 'F%RMU     ',SIZE(f%RMU ),shape(f%RMU ) 
ALLOCATE(f%RW(r%NDGL))
IF (llp2) WRITE(nout,9) 'F%RW      ',SIZE(f%RW  ),shape(f%RW  ) 
ALLOCATE(f%R1MU2(r%NDGL))
IF (llp2) WRITE(nout,9) 'F%R1MU2   ',SIZE(f%R1MU2),shape(f%R1MU2 ) 
ALLOCATE(f%RACTHE(r%NDGL))
IF (llp2) WRITE(nout,9) 'F%RACTHE  ',SIZE(f%RACTHE),shape(f%RACTHE ) 

!CALL GSTATS(1801,0) ! MPL 4.12.08
DO jnm=1,d%NSPOLEGL
  f%RPNM(r%NLEI3,jnm) = 0.0_jprb
ENDDO
!CALL GSTATS(1801,1)     ! MPL 4.12.08

!     ------------------------------------------------------------------

!*       3.1   Gaussian latitudes and weights
!CALL SUGAW(R%NDGL,F%RMU,DLRMU,F%RW)

!*       3.2   Computes related arrays

DO jgl=1,r%NDGL
  f%R1MU2(jgl) = REAL(1.0_JPRB-DLRMU(JGL)*DLRMU(JGL),JPKS)
  f%RACTHE(jgl) = REAL(1.0_jprb/sqrt(1.0_jprb-dlrmu(jgl)*dlrmu(jgl))/&
   &REAL(RA,JPKD),JPKS)
ENDDO

!*       3.3   Working arrays
DO jn=3,r%NTMAX+1
  DO jm=2,jn-1
    dlc(jm,jn) = dc(jn,jm)
    dld(jm,jn) = dd(jn,jm)
    dle(jm,jn) = de(jn,jm)
  ENDDO
ENDDO

DO jn=1,r%NTMAX+1
  dla(jn) = sqrt(REAL(2*jn+1,jpkd))
  dlb(jn) = sqrt(REAL(2*jn+1,jpkd)/REAL(JN*(JN+1),JPKD))
  dlf(jn) = REAL(2*jn-1,jpkd)/REAL(jn,jpkd)
  dlg(jn) = REAL(jn-1,jpkd)/REAL(jn,jpkd)
  dlh(jn) = sqrt(REAL(2*jn+1,jpkd)/REAL(2*jn,jpkd))
  dli(jn) = REAL(jn,jpkd)
ENDDO

!CALL GSTATS(1801,0)  ! MPL 4.12.08
DO jgl=d%NLATLS(mysetw),d%NLATLE(mysetw)
  dlpol(:,:) = 0.0_jprb
  CALL supol(r%NTMAX+1,dlrmu(jgl),dlpol,dla,dlb,dlc,dld,dle,dlf,dlg,dlh,dli)
  inm = 0
  iglloc = jgl - d%NLATLS(mysetw) + 1
  DO jm=0,r%NSMAX
    DO jn=r%NTMAX+1,jm,-1
      inm = inm+1
      zpnmg(inm,iglloc) = REAL(DLPOL(JM,JN),JPKS)
    ENDDO
  ENDDO
ENDDO
!CALL GSTATS(1801,1) !  MPL 4.12.08
!CALL GSTATS(140,1)  ! MPL 4.12.08

!CALL GSTATS(190,0)  ! MPL 4.12.08
CALL sutrle(zpnmg)
!CALL GSTATS(190,1)  ! MPL 4.12.08

icount = 0
DO jmloc=1,d%NUMP
  im = d%MYMS(jmloc)
  DO jn=im,r%NTMAX+2
    icount = icount+1
  ENDDO
ENDDO

ALLOCATE(f%REPSNM(icount))
IF (llp2) WRITE(nout,9) 'F%REPSNM  ',SIZE(f%REPSNM ),shape(f%REPSNM ) 

icount = 0
DO jmloc=1,d%NUMP
  im = d%MYMS(jmloc)
  DO jn=im,r%NTMAX+2
    icount = icount+1
    f%REPSNM(icount) = REAL(SQRT(REAL(JN*JN-IM*IM,JPKD)/&
     &REAL(4*JN*JN-1,JPKD)),JPKS)
  ENDDO
ENDDO

ALLOCATE(f%RN(-1:r%NTMAX+3))
IF (llp2) WRITE(nout,9) 'F%RN      ',SIZE(f%RN     ),shape(f%RN     ) 
ALLOCATE(f%RLAPIN(-1:r%NSMAX+2))
IF (llp2) WRITE(nout,9) 'F%RLAPIN  ',SIZE(f%RLAPIN ),shape(f%RLAPIN ) 
ALLOCATE(f%NLTN(-1:r%NTMAX+3))
IF (llp2) WRITE(nout,9) 'F%NLTN    ',SIZE(f%NLTN ),shape(f%NLTN ) 

DO jn=-1,r%NTMAX+3
  f%RN(jn) = REAL(jn,jprb)
  f%NLTN(jn) = r%NTMAX+2-jn
ENDDO
f%RLAPIN(:) = 0.0_jprb
f%RLAPIN(0) = 0._jprb
f%RLAPIN(-1) = 0.0_jprb
DO jn=1,r%NSMAX+2
  f%RLAPIN(jn)=REAL(-(REAL(RA,JPKD)*REAL(RA,JPKD))/REAL(JN*(JN+1),JPKD),JPKS)
ENDDO

DEALLOCATE(zpnmg)

!     ------------------------------------------------------------------
9 FORMAT(1x,'ARRAY ',a10,' ALLOCATED ',8i8)

END SUBROUTINE suleg
END MODULE suleg_mod