source: LMDZ5/trunk/libf/phylmd/rrtm/suleg_mod.F90 @ 2930

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