source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/rrtm/suleg_mod.F90 @ 3331

Last change on this file since 3331 was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 7.6 KB
Line 
1MODULE SULEG_MOD
2CONTAINS
3SUBROUTINE SULEG
4
5USE PARKIND1  ,ONLY : JPIM     ,JPRB
6USE PARKIND2  ,ONLY : JPRH
7
8USE TPM_GEN
9USE TPM_DIM
10USE TPM_CONSTANTS
11USE TPM_DISTR
12USE TPM_FIELDS
13
14!USE SUGAW_MOD
15USE SUPOL_MOD
16USE SUTRLE_MOD
17
18#ifdef DOC
19
20!**** *SULEG * - initialize the Legendre polynomials
21
22!     Purpose.
23!     --------
24!           Initialize COMMON YOMLEG
25
26!**   Interface.
27!     ----------
28!        *CALL* *SULEG*
29
30!        Explicit arguments :
31!        --------------------
32
33!        Implicit arguments :
34!        --------------------
35!              COMMON YOMLEG
36
37!     Method.
38!     -------
39!        See documentation
40
41!     Externals.
42!     ----------
43!        SUGAW  (Gaussian latitudes)
44!        SUPOLM (polynomials)
45!        LFI routines for external IO's
46!        Called by SUGEM.
47
48!     Reference.
49!     ----------
50!        ECMWF Research Department documentation of the IFS
51
52!     Author.
53!     -------
54!        Mats Hamrud and Philippe Courtier  *ECMWF*
55
56!     Modifications.
57!     --------------
58!        Original : 87-10-15
59!        MODIFICATION : 91-04  J.M. Piriou:
60!                       - Read gaussian latitudes and PNM on LFI
61!                       - If file missing, computes
62!                       91-04 M.Hamrud:
63!                       - IO Scheme introduced
64!        MODIFICATION : 91-07-03 P.Courtier suppress derivatives
65!        MODIFICATION : 91-07-03 P.Courtier computes RATATH and RACTHE
66!        MODIFICATION : 91-07-03 P.Courtier change upper limit (NSMAX+1)
67!        MODIFICATION : 91-07-03 P.Courtier change ordering
68!     Order of the PNM in the file, as in the model :
69!         - increasing wave numbers m
70!         - for a given m, from n=NSMAX+1 to m
71!        MODIFICATION : 92-07-02 R. Bubnova: shift RATATH calculation
72!                                            to SUGEM1
73!        MODIFICATION : 92-12-17 P.Courtier multitask computations
74!        Modified by R. EL Khatib : 93-04-02 Set-up defaults controled by LECMWF
75!        MODIFICATION : 93-03-19 D.Giard : n <= NTMAX
76!        K. YESSAD    : 93-05-11 : DLMU --> global array DRMU(NDGSA:NDGEN).
77!                       (not stored currently on LFI files).
78!        MODIFICATION : 94-02-03 R. El Khatib : subroutine SULEG2 to write out
79!                       the Leg. polynomials on workfile or LFI file
80!        Modification : 94-08-31 M. Tolstykh: Setup for CUD interpolation
81!        Modified by K. YESSAD (MARCH 1995): Extra-latitudes computations
82!                 according to value of NDGSUR and LRPOLE only.
83!                 + change fancy loop numbering.
84!        Modified 98-08-10 by K. YESSAD: removal of LRPOLE option.
85!          - removal of LRPOLE in YOMCT0.
86!          - removal of code under LRPOLE.
87!     ------------------------------------------------------------------
88
89#endif
90
91
92IMPLICIT NONE
93
94
95!     ------------------------------------------------------------------
96REAL(KIND=JPRB),ALLOCATABLE :: ZPNMG(:,:)
97
98REAL(KIND=JPRH) :: DLRMU(R%NDGL)
99REAL(KIND=JPRH) :: DLC(0:R%NTMAX+1,0:R%NTMAX+1)
100REAL(KIND=JPRH) :: DLD(0:R%NTMAX+1,0:R%NTMAX+1)
101REAL(KIND=JPRH) :: DLE(0:R%NTMAX+1,0:R%NTMAX+1)
102REAL(KIND=JPRH) :: DLA(0:R%NTMAX+1),DLB(0:R%NTMAX+1),DLF(0:R%NTMAX+1)
103REAL(KIND=JPRH) :: DLG(0:R%NTMAX+1),DLH(0:R%NTMAX+1),DLI(0:R%NTMAX+1)
104REAL(KIND=JPRH) :: DLPOL(0:R%NTMAX+1,0:R%NTMAX+1)
105!     ------------------------------------------------------------------
106
107INTEGER(KIND=JPIM), PARAMETER :: JPKS=KIND(ZPNMG)
108INTEGER(KIND=JPIM), PARAMETER :: JPKD=KIND(DLG)
109
110!     ------------------------------------------------------------------
111REAL(KIND=JPRH) :: DA,DC,DD,DE
112INTEGER(KIND=JPIM) :: KKN, KKM
113
114!     LOCAL
115INTEGER(KIND=JPIM) :: IGLLOC, INM, IM , ICOUNT,&
116             &JGL,  JM, JMLOC, JN, JNM
117
118
119LOGICAL :: LLP1,LLP2
120
121
122DC(KKN,KKM)=SQRT( (REAL(2*KKN+1,JPKD)*REAL(KKN+KKM-1,JPKD)&
123                   &*REAL(KKN+KKM-3,JPKD))&
124                &/ (REAL(2*KKN-3,JPKD)*REAL(KKN+KKM,JPKD)&
125                   &*REAL(KKN+KKM-2,JPKD)) )
126DD(KKN,KKM)=SQRT( (REAL(2*KKN+1,JPKD)*REAL(KKN+KKM-1,JPKD)&
127                   &*REAL(KKN-KKM+1,JPKD))&
128                &/ (REAL(2*KKN-1,JPKD)*REAL(KKN+KKM,JPKD)&
129                  &*REAL(KKN+KKM-2,JPKD)) )
130DE(KKN,KKM)=SQRT( (REAL(2*KKN+1,JPKD)*REAL(KKN-KKM,JPKD))&
131                &/ (REAL(2*KKN-1,JPKD)*REAL(KKN+KKM,JPKD)) )
132DA(KKN,KKM)=SQRT( (REAL(2*KKN+1,JPKD)*REAL(KKN-KKM,JPKD)&
133                   &*REAL(KKN+KKM,JPKD))&
134                &/  REAL(2*KKN-1,JPKD) )
135
136!     ------------------------------------------------------------------
137ALLOCATE(ZPNMG(R%NSPOLEG,D%NLEI3D))
138
139!*       0.    Some initializations.
140!              ---------------------
141
142LLP1 = NPRINTLEV>0
143LLP2 = NPRINTLEV>1
144IF(LLP1) WRITE(NOUT,*) '=== ENTER ROUTINE SULEG ==='
145
146!CALL GSTATS(140,0)  !MPL 4.12.08
147ALLOCATE(F%RPNM(R%NLEI3,D%NSPOLEGL))
148IF (LLP2) WRITE(NOUT,9) 'F%RPNM    ',SIZE(F%RPNM),SHAPE(F%RPNM)
149ALLOCATE(F%RMU(R%NDGL))
150IF (LLP2) WRITE(NOUT,9) 'F%RMU     ',SIZE(F%RMU ),SHAPE(F%RMU )
151ALLOCATE(F%RW(R%NDGL))
152IF (LLP2) WRITE(NOUT,9) 'F%RW      ',SIZE(F%RW  ),SHAPE(F%RW  )
153ALLOCATE(F%R1MU2(R%NDGL))
154IF (LLP2) WRITE(NOUT,9) 'F%R1MU2   ',SIZE(F%R1MU2),SHAPE(F%R1MU2 )
155ALLOCATE(F%RACTHE(R%NDGL))
156IF (LLP2) WRITE(NOUT,9) 'F%RACTHE  ',SIZE(F%RACTHE),SHAPE(F%RACTHE )
157
158!CALL GSTATS(1801,0) ! MPL 4.12.08
159DO JNM=1,D%NSPOLEGL
160  F%RPNM(R%NLEI3,JNM) = 0.0_JPRB
161ENDDO
162!CALL GSTATS(1801,1)     ! MPL 4.12.08
163
164!     ------------------------------------------------------------------
165
166!*       3.1   Gaussian latitudes and weights
167!CALL SUGAW(R%NDGL,F%RMU,DLRMU,F%RW)
168
169!*       3.2   Computes related arrays
170
171DO JGL=1,R%NDGL
172  F%R1MU2(JGL) = REAL(1.0_JPRB-DLRMU(JGL)*DLRMU(JGL),JPKS)
173  F%RACTHE(JGL) = REAL(1.0_JPRB/SQRT(1.0_JPRB-DLRMU(JGL)*DLRMU(JGL))/&
174   &REAL(RA,JPKD),JPKS)
175ENDDO
176
177!*       3.3   Working arrays
178DO JN=3,R%NTMAX+1
179  DO JM=2,JN-1
180    DLC(JM,JN) = DC(JN,JM)
181    DLD(JM,JN) = DD(JN,JM)
182    DLE(JM,JN) = DE(JN,JM)
183  ENDDO
184ENDDO
185
186DO JN=1,R%NTMAX+1
187  DLA(JN) = SQRT(REAL(2*JN+1,JPKD))
188  DLB(JN) = SQRT(REAL(2*JN+1,JPKD)/REAL(JN*(JN+1),JPKD))
189  DLF(JN) = REAL(2*JN-1,JPKD)/REAL(JN,JPKD)
190  DLG(JN) = REAL(JN-1,JPKD)/REAL(JN,JPKD)
191  DLH(JN) = SQRT(REAL(2*JN+1,JPKD)/REAL(2*JN,JPKD))
192  DLI(JN) = REAL(JN,JPKD)
193ENDDO
194
195!CALL GSTATS(1801,0)  ! MPL 4.12.08
196DO JGL=D%NLATLS(MYSETW),D%NLATLE(MYSETW)
197  DLPOL(:,:) = 0.0_JPRB
198  CALL SUPOL(R%NTMAX+1,DLRMU(JGL),DLPOL,DLA,DLB,DLC,DLD,DLE,DLF,DLG,DLH,DLI)
199  INM = 0
200  IGLLOC = JGL - D%NLATLS(MYSETW) + 1
201  DO JM=0,R%NSMAX
202    DO JN=R%NTMAX+1,JM,-1
203      INM = INM+1
204      ZPNMG(INM,IGLLOC) = REAL(DLPOL(JM,JN),JPKS)
205    ENDDO
206  ENDDO
207ENDDO
208!CALL GSTATS(1801,1) !  MPL 4.12.08
209!CALL GSTATS(140,1)  ! MPL 4.12.08
210
211!CALL GSTATS(190,0)  ! MPL 4.12.08
212CALL SUTRLE(ZPNMG)
213!CALL GSTATS(190,1)  ! MPL 4.12.08
214
215ICOUNT = 0
216DO JMLOC=1,D%NUMP
217  IM = D%MYMS(JMLOC)
218  DO JN=IM,R%NTMAX+2
219    ICOUNT = ICOUNT+1
220  ENDDO
221ENDDO
222
223ALLOCATE(F%REPSNM(ICOUNT))
224IF (LLP2) WRITE(NOUT,9) 'F%REPSNM  ',SIZE(F%REPSNM ),SHAPE(F%REPSNM )
225
226ICOUNT = 0
227DO JMLOC=1,D%NUMP
228  IM = D%MYMS(JMLOC)
229  DO JN=IM,R%NTMAX+2
230    ICOUNT = ICOUNT+1
231    F%REPSNM(ICOUNT) = REAL(SQRT(REAL(JN*JN-IM*IM,JPKD)/&
232     &REAL(4*JN*JN-1,JPKD)),JPKS)
233  ENDDO
234ENDDO
235
236ALLOCATE(F%RN(-1:R%NTMAX+3))
237IF (LLP2) WRITE(NOUT,9) 'F%RN      ',SIZE(F%RN     ),SHAPE(F%RN     )
238ALLOCATE(F%RLAPIN(-1:R%NSMAX+2))
239IF (LLP2) WRITE(NOUT,9) 'F%RLAPIN  ',SIZE(F%RLAPIN ),SHAPE(F%RLAPIN )
240ALLOCATE(F%NLTN(-1:R%NTMAX+3))
241IF (LLP2) WRITE(NOUT,9) 'F%NLTN    ',SIZE(F%NLTN ),SHAPE(F%NLTN )
242
243DO JN=-1,R%NTMAX+3
244  F%RN(JN) = REAL(JN,JPRB)
245  F%NLTN(JN) = R%NTMAX+2-JN
246ENDDO
247F%RLAPIN(:) = 0.0_JPRB
248F%RLAPIN(0) = 0._JPRB
249F%RLAPIN(-1) = 0.0_JPRB
250DO JN=1,R%NSMAX+2
251  F%RLAPIN(JN)=REAL(-(REAL(RA,JPKD)*REAL(RA,JPKD))/REAL(JN*(JN+1),JPKD),JPKS)
252ENDDO
253
254DEALLOCATE(ZPNMG)
255
256!     ------------------------------------------------------------------
2579 FORMAT(1X,'ARRAY ',A10,' ALLOCATED ',8I8)
258
259END SUBROUTINE SULEG
260END MODULE SULEG_MOD
Note: See TracBrowser for help on using the repository browser.