source: LMDZ5/branches/testing/libf/phymar/radaca.F90 @ 3946

Last change on this file since 3946 was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
File size: 18.5 KB
RevLine 
[2089]1!OPTIONS XOPT(HSFUN)
2SUBROUTINE RADACA &
3 &( KIDIA , KFDIA , KLON , KTDIA , KLEV &
4 &, PAPRS , PGELAM, PSIN  , PCLON, PSLON , PTH &
5 &, PAER  , POZON &
6 &)
7
8!***********************************************************************
9! CAUTION: THIS ROUTINE WORKS ONLY ON A NON-ROTATED, UNSTRETCHED GRID
10!***********************************************************************
11
12!**** *RADACA  - COMPUTES DISTRIBUTION OF AEROSOLS AND OZONE
13
14!     PURPOSE.
15!     --------
16
17!**   INTERFACE.
18!     ----------
19!        CALL *RADACA* FROM *RADINT*
20
21!        EXPLICIT ARGUMENTS :
22!        --------------------
23!     ==== INPUTS ===
24!     ==== OUTPUTS ===
25
26!        IMPLICIT ARGUMENTS :   NONE
27!        --------------------
28
29!     METHOD.
30!     -------
31
32
33!     EXTERNALS.
34!     ----------
35
36!          NONE
37
38!     REFERENCE.
39!     ----------
40
41!        SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND
42!        ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S"
43
44!     AUTHOR.
45!     -------
46!     J.-J. MORCRETTE  E.C.M.W.F.    91/03/15
47
48!     MODIFICATIONS.
49!     --------------
50!     J.-J. MORCRETTE  E.C.M.W.F.    93/03/15   OPERATIONAL CLIMATOLOGY
51!     JJMorcrette 98-12-21 GISS volcanic aerosol climatology
52!     JJMorcrette 99-09    monthly climatology of tropospheric aerosols
53
54!     ADAPTATION TO zahir.idris.fr
55!     ----------------------------
56!     HGallée     04-11-17 Definition of the arguments 2 & 3 of LEGTRI
57!-----------------------------------------------------------------------
58
59#include "tsmbkind.h"
60
61USE YOMCST   , ONLY : R        ,RPI
62USE YOEAERD  , ONLY : CVDAES   ,CVDAEL   ,CVDAEU   ,CVDAED   ,&
63            &RCAEOPS  ,RCAEOPL  ,RCAEOPU  ,RCAEOPD  ,RCTRBGA  ,&
64            &RCVOBGA  ,RCSTBGA  ,RCTRPT   ,RAESC    ,RAESS    ,&
65            &RAELC    ,RAELS    ,RAEUC    ,RAEUS    ,RAEDC    ,&
66            &RAEDS
67USE YOEOZOC  , ONLY : COZQC    ,COZQS    ,COZHC    ,COZHS
68USE YOERAD   , ONLY : LHVOLCA  ,LNEWAER
69USE YOEAERC  , ONLY : RSINCT   ,RSINCV   ,REPAER   ,&
70            &RTAEBC  ,RTAEOR   ,RTAESD   ,RTAESS   ,RTAESU   , &
71            &RTAEVO
72
73
74IMPLICIT NONE
75
76
77!     DUMMY INTEGER SCALARS
78INTEGER_M :: KFDIA
79INTEGER_M :: KIDIA
80INTEGER_M :: KLEV
81INTEGER_M :: KLON
82INTEGER_M :: KTDIA
83
84INTEGER_M :: KCF, KRINT, KSHIFT
85
86INTEGER_M :: KCP__RADACA
87INTEGER_M :: KDIM_RADACA
88
89!     -----------------------------------------------------------------
90
91!*       0.1   ARGUMENTS.
92!              ----------
93
94REAL_B :: PAPRS(KLON,KLEV+1), PGELAM(KLON), PSIN(KLON)&
95 &,   PCLON(KLON)       , PSLON(KLON) , PTH(KLON,KLEV+1)
96
97REAL_B :: PAER(KLON,6,KLEV),POZON(KLON,KLEV)
98!     -----------------------------------------------------------------
99
100!*       0.2   LOCAL ARRAYS.
101!              -------------
102
103INTEGER_M :: IINLA1(KLON), IINLA2(KLON)
104INTEGER_M :: IINLO1(KLON), IINLO2(KLON)
105
106REAL_B :: ZAED  (KLON), ZAEL  (KLON), ZAES  (KLON), ZAEU  (KLON)
107REAL_B :: ZAEQDN(KLON), ZAEQDO(KLON), ZAEQLN(KLON), ZAEQLO(KLON)
108REAL_B :: ZAEQSN(KLON), ZAEQSO(KLON), ZAEQUN(KLON), ZAEQUO(KLON)
109
110REAL_B :: ZAERBC(KLON), ZAEROR(KLON), ZAERSD(KLON)
111REAL_B :: ZAERSS(KLON), ZAERSU(KLON), ZAERVO(KLON)
112
113REAL_B :: ZAETRN(KLON),ZAETRO(KLON)
114
115REAL_B :: ZALP(66)
116REAL_B :: ZDPN(KLON)  , ZDPO(KLON)
117REAL_B :: ZFAED(21)    , ZFAEL(21)    , ZFAES(21)    , ZFAEU(21)
118REAL_B :: ZFOZQ(11)    , ZFOZH(11)
119REAL_B :: ZGRTH(KLON)
120REAL_B :: ZLON(KLON)  , ZLONR(72)    , ZNLO1(KLON) , ZNLO2(KLON)
121
122REAL_B :: ZOZH (KLON) , ZOZQ (KLON)
123REAL_B :: ZQOFN(KLON) , ZQOFO(KLON)
124REAL_B :: ZSILAT(KLON), ZSINR(46)
125
126!     LOCAL INTEGER SCALARS
127INTEGER_M :: IL, IMM, IMNC, IMNS, INLA, INLA1, INLA2, INLO1, INLO2, &
128          &ITOTPT, JK, JL, JLR, JMM, &
129          &JNN, NLATR, NLONR, JAER, JEND, JIL, JJL, IPRINT, ITOT
130
131!     LOCAL REAL SCALARS
132REAL_B :: ZAETR, ZCOS1, ZCOS10, ZCOS2, ZCOS3, ZCOS4,&
133          &ZCOS5, ZCOS6, ZCOS7, ZCOS8, ZCOS9, ZCPHN3, &
134          &ZCPHO3, ZDPNMO, ZGRIDR, ZLATR, ZSDPN3, ZSDPO3, &
135          &ZSIN, ZSIN1, ZSIN10, ZSIN2, ZSIN3, ZSIN4, &
136          &ZSIN5, ZSIN6, ZSIN7, ZSIN8, ZSIN9
137REAL_B :: ZAERBC1, ZAERBC2, ZAEROR1, ZAEROR2, ZAERSD1, ZAERSD2, &
138          &ZAERSS1, ZAERSS2, ZAERSU1, ZAERSU2
139REAL_B :: ZDLONR
140
141!     ------------------------------------------------------------------
142
143!*         1.     "NEW AEROSOL DISTRIBUTION" PARAMETERS COMPUTATIONS
144!                 --------------------------------------------------
145
146
147!*         1.1    VOLCANIC AEROSOL DISTRIBUTION PARAMETERS
148!                 ----------------------------------------
149
150!                 GISS CLIMATOLOGY
151!                 ----------------
152
153KCF=0
154KSHIFT=0
155KRINT=1
156
157IF (LHVOLCA) THEN
158  NLATR=46
159  ZGRIDR=180._JPRB/(NLATR-1)
160  DO JLR=1,NLATR
161    ZLATR=90._JPRB-(JLR-1)*ZGRIDR
162    ZSINR(JLR)=SIN(ZLATR*RPI/180._JPRB)
163  ENDDO
164
165  IL=KSHIFT
166  DO JL=KIDIA,KFDIA
167    IL=IL+1
168    INLA=0
169    ZSILAT(IL)=-9999._JPRB
170    ZSIN=PSIN(JL)
171    DO JLR=1,NLATR-1
172      IF (ZSIN <= ZSINR(JLR) .AND. ZSIN > ZSINR(JLR+1)) THEN
173        INLA=JLR
174        ZSILAT(IL)=(ZSIN-ZSINR(INLA))/(ZSINR(INLA+1)-ZSINR(INLA))
175        ZAERVO(IL)=RTAEVO(INLA)+ZSILAT(IL)*(RTAEVO(INLA+1)-RTAEVO(INLA))
176      ENDIF
177    ENDDO
178    IF (ZSIN <= ZSINR(NLATR-1) .AND. ZSIN >= ZSINR(NLATR))THEN
179      INLA=NLATR
180      ZSILAT(IL)=(ZSIN-ZSINR(INLA-1))/(ZSINR(INLA)-ZSINR(INLA-1))
181      ZAERVO(IL)=RTAEVO(INLA-1)&
182       &+ZSILAT(IL)*(RTAEVO(INLA)-RTAEVO(INLA-1))
183    ENDIF
184    IF (INLA == 0) THEN
185!      CALL ABOR1(' Problem with lat. interpolation in radaca!')
186      STOP ' Problem with lat. interpolation in radaca!'
187    ENDIF
188  ENDDO
189
190!                 TANRE ET AL. CLIMATOLOGY
191!                 ------------------------
192
193ELSE
194  IL = KSHIFT
195  DO JL=KIDIA,KFDIA
196    IL = IL+1
197    ZAERVO(IL)=RCVOBGA
198  ENDDO
199ENDIF
200ITOTPT=IL
201
202
203
204!*         1.2    TROPOSPHERIC AEROSOL DISTRIBUTION PARAMETERS
205!                 --------------------------------------------
206
207IF (LNEWAER) THEN
208!  print *,'LNEWAER= ',LNEWAER
209!  print *,' inputs        SINLAT            LONGITUDE'
210!  DO JL=KIDIA,KFDIA
211!    print 9001,JL,PSIN(JL),PGELAM(JL)
2129001    format(1x,'RADACA ',1x,I5,1x,2E15.8)
213!  END DO         
214
215!-- latitude       
216  NLATR=46
217  ZGRIDR=180._JPRB/(NLATR-1)
218  DO JLR=1,NLATR
219    ZLATR=90._JPRB-(JLR-1)*ZGRIDR
220    ZSINR(JLR)=SIN(ZLATR*RPI/180._JPRB)
221  END DO
222  NLONR=72
223  ZDLONR=2._JPRB*RPI/NLONR
224  DO JLR=1,NLONR
225    ZLONR(JLR)=(JLR-1)*ZDLONR
226  END DO
227
228!  print *,'Reference grid for Tegen climatology'       
229!  print 9121,(ZSINR(JLR),JLR=1,NLATR)
2309121     format(1x,'ZSINR ',8E15.7)
231!  print 9122,(ZLONR(JLR),JLR=1,NLONR)
2329122     format(1x,'ZLONR ',8E15.7)
233           
234  IL=KSHIFT
235  DO JL=KIDIA,KFDIA         
236    IL=IL+1
237    IINLA1(IL)=0
238    IINLA2(IL)=0
239    ZSILAT(IL)=-9999._JPRB
240    ZSIN=PSIN(JL)
241    DO JLR=1,NLATR-1
242      IF (ZSIN <= ZSINR(JLR) .AND. ZSIN > ZSINR(JLR+1)) THEN
243        INLA=JLR
244        IINLA1(IL)=JLR
245        IINLA2(IL)=JLR+1
246        ZSILAT(IL)=(ZSIN-ZSINR(INLA))/(ZSINR(INLA+1)-ZSINR(INLA))
247      ENDIF
248    ENDDO
249    IF (ZSIN <= ZSINR(NLATR-1) .AND. ZSIN >= ZSINR(NLATR))THEN
250      INLA=NLATR
251      IINLA1(IL)=NLATR-1
252      IINLA2(IL)=NLATR
253      ZSILAT(IL)=(ZSIN-ZSINR(INLA-1))/(ZSINR(INLA)-ZSINR(INLA-1))
254    END IF   
255    IF (INLA.EQ.0) THEN
256!      CALL ABOR1(' Problem with lat. interpolation in radaca!')
257      STOP ' Problem with lat. interpolation in radaca!'
258    ENDIF
259!    print 9123,JL,IL,PSIN(JL),INLA,ZSINR(INLA),ZSILAT(IL)
2609123     format(1x,'Interp.Latit.',2I4,F10.7,I4,2F10.7)         
261  END DO 
262
263!-- longitude
264  IL=KSHIFT
265  DO JL=KIDIA,KFDIA
266    IL=IL+1
267    IINLO1(IL)=0
268    IINLO2(IL)=0
269    ZLON(IL)=-9999.
270    DO JLR=1,71
271      IF (PGELAM(JL) < ZLONR(JLR+1) .AND. PGELAM(JL) >= ZLONR(JLR)) &
272     &  THEN
273        IINLO1(IL)=JLR
274        IINLO2(IL)=JLR+1
275        ZNLO1(IL)=ZLONR(JLR)
276        ZNLO2(IL)=ZLONR(JLR+1)
277      END IF
278    END DO     
279    IF (PGELAM(JL) >= ZLONR(72)) THEN
280      IINLO1(IL)=72
281      IINLO2(IL)= 1
282      ZNLO1(IL)=ZLONR(72)
283!      ZNLO2(IL)=ZLONR(72)+2.*RPI
284      ZNLO2(IL)=ZLONR(72)+ZDLONR
285    ENDIF
286! Martin control
287!    print 9124,JL,IL,PGELAM(JL),IINLO1(IL),IINLO2(IL),ZNLO1(IL),ZNLO2(IL)
288!9124     format(1x,'Interp.Longit0.',2I4,F10.7,2I5,2F10.7)         
289  END DO 
290 
291  IL=KSHIFT
292  DO JL=KIDIA,KFDIA       
293    IL=IL+1
294    IF (IINLO1(IL).EQ.0 .OR. IINLO2(IL).EQ.0) THEN
295!      CALL ABOR1(' Problem with long. interpolation in radaca!')
296      STOP ' Problem with long. interpolation in radaca!'
297    ENDIF
298    ZLON(IL)=(PGELAM(JL)-ZNLO1(IL))/(ZNLO2(IL)-ZNLO1(IL))
299    INLO1=IINLO1(IL)
300    INLO2=IINLO2(IL)
301    INLA1=IINLA1(IL)
302    INLA2=IINLA2(IL)
303   
304    ZAERBC1=RTAEBC(INLO1,INLA1) &
305     &      +ZSILAT(IL)*(RTAEBC(INLO1,INLA2)-RTAEBC(INLO1,INLA1)) 
306    ZAERBC2=RTAEBC(INLO2,INLA1) &
307     &      +ZSILAT(IL)*(RTAEBC(INLO2,INLA2)-RTAEBC(INLO2,INLA1)) 
308    ZAERBC(IL)=ZAERBC1+ZLON(IL)*(ZAERBC2-ZAERBC1)
309       
310    ZAEROR1=RTAEOR(INLO1,INLA1) &
311     &      +ZSILAT(IL)*(RTAEOR(INLO1,INLA2)-RTAEOR(INLO1,INLA1)) 
312    ZAEROR2=RTAEOR(INLO2,INLA1) &
313     &      +ZSILAT(IL)*(RTAEOR(INLO2,INLA2)-RTAEOR(INLO2,INLA1)) 
314    ZAEROR(IL)=ZAEROR1+ZLON(IL)*(ZAEROR2-ZAEROR1)
315       
316    ZAERSD1=RTAESD(INLO1,INLA1) &
317     &      +ZSILAT(IL)*(RTAESD(INLO1,INLA2)-RTAESD(INLO1,INLA1)) 
318    ZAERSD2=RTAESD(INLO2,INLA1) &
319     &      +ZSILAT(IL)*(RTAESD(INLO2,INLA2)-RTAESD(INLO2,INLA1)) 
320    ZAERSD(IL)=ZAERSD1+ZLON(IL)*(ZAERSD2-ZAERSD1)
321       
322    ZAERSS1=RTAESS(INLO1,INLA1) &
323     &      +ZSILAT(IL)*(RTAESS(INLO1,INLA2)-RTAESS(INLO1,INLA1)) 
324    ZAERSS2=RTAESS(INLO2,INLA1) &
325     &      +ZSILAT(IL)*(RTAESS(INLO2,INLA2)-RTAESS(INLO2,INLA1)) 
326    ZAERSS(IL)=ZAERSS1+ZLON(IL)*(ZAERSS2-ZAERSS1)
327       
328    ZAERSU1=RTAESU(INLO1,INLA1) &
329     &      +ZSILAT(IL)*(RTAESU(INLO1,INLA2)-RTAESU(INLO1,INLA1)) 
330    ZAERSU2=RTAESU(INLO2,INLA1) &
331     &      +ZSILAT(IL)*(RTAESU(INLO2,INLA2)-RTAESU(INLO2,INLA1)) 
332    ZAERSU(IL)=ZAERSU1+ZLON(IL)*(ZAERSU2-ZAERSU1)
333         
334!    print 9125,JL,IL,PSIN(JL),PGELAM(JL),ZSILAT(IL) &
335!     &      ,RTAESU(INLO1,INLA2),RTAESU(INLO1,INLA1),ZAERSU1 &
336!     &      ,RTAESU(INLO2,INLA2),RTAESU(INLO2,INLA1),ZAERSU2 &
337!     &                    ,INLA1,INLO1,INLO2,INLA2
338!9125     format(1x,'Interp.Longit1.',2I4,9F10.7,4I5)         
339!    print 9126,JL,IL,PSIN(JL),PGELAM(JL),ZSILAT(IL),ZLON(IL) &
340!     &      ,ZNLO1(IL),ZNLO2(IL),INLA1,INLO1,INLO2,INLA2
341!9126     format(1x,'Interp.Longit2.',2I4,6F10.7,4I5)         
342!    print 9127,JL,IL,ZAERBC(IL),ZAEROR(IL),ZAERSD(IL),ZAERSS(IL) &
343!     &      ,ZAERSU(IL)
344!9127     format(1x,'Interp.Longit3.',2I4,5F10.7)         
345  END DO       
346END IF
347
348!     ------------------------------------------------------------------   
349
350!*       2.     OZONE
351!               -----
352
353ZSIN=PSIN(KIDIA)
354
355!*       2.1     CALL TO LEGTRI.
356!                ---------------
357KCP__RADACA= 6
358KDIM_RADACA=66
359!***
360CALL LEGTRI (ZSIN,KCP__RADACA,KDIM_RADACA,ZALP)
361!***
362
363!*       2.2     LEGENDRE TRANSFORM FOR OZONE.
364!                -----------------------------
365
366DO JMM=1,11
367  ZFOZQ(JMM)=_ZERO_
368  ZFOZH(JMM)=_ZERO_
369ENDDO
370IMM=0
371IMNC=0
372IMNS=0
373DO JMM=1,6
374  IMM=IMM+1
375  DO JNN=JMM,6
376    IMNC=IMNC+1
377    ZFOZQ(IMM)=ZFOZQ(IMM)+ZALP(IMNC)*COZQC(IMNC)
378    ZFOZH(IMM)=ZFOZH(IMM)+ZALP(IMNC)*COZHC(IMNC)
379  ENDDO
380  IF(JMM /= 1) THEN
381    IMM=IMM+1
382    DO JNN=JMM,6
383      IMNS=IMNS+1
384      ZFOZQ(IMM)=ZFOZQ(IMM)+ZALP(IMNS+6)*COZQS(IMNS)
385      ZFOZH(IMM)=ZFOZH(IMM)+ZALP(IMNS+6)*COZHS(IMNS)
386    ENDDO
387  ENDIF
388ENDDO
389
390!*       2.3     FOURIER TRANSFORM FOR OZONE.
391!                ----------------------------
392
393IL=KSHIFT
394DO JL=KIDIA,KFDIA
395  IL=IL+1
396  ZCOS1=PCLON(JL)
397  ZSIN1=PSLON(JL)
398  ZCOS2=ZCOS1*ZCOS1-ZSIN1*ZSIN1
399  ZSIN2=ZSIN1*ZCOS1+ZCOS1*ZSIN1
400  ZCOS3=ZCOS2*ZCOS1-ZSIN2*ZSIN1
401  ZSIN3=ZSIN2*ZCOS1+ZCOS2*ZSIN1
402  ZCOS4=ZCOS3*ZCOS1-ZSIN3*ZSIN1
403  ZSIN4=ZSIN3*ZCOS1+ZCOS3*ZSIN1
404  ZCOS5=ZCOS4*ZCOS1-ZSIN4*ZSIN1
405  ZSIN5=ZSIN4*ZCOS1+ZCOS4*ZSIN1
406  ZOZQ(IL)=&
407   &ZFOZQ(1)+_TWO_*(ZFOZQ(2)*ZCOS1+ZFOZQ(3)*ZSIN1+ZFOZQ(4)*ZCOS2 &
408   &+ZFOZQ(5)*ZSIN2+ZFOZQ(6)*ZCOS3+ZFOZQ(7)*ZSIN3+ZFOZQ(8)&
409   &*ZCOS4+ZFOZQ(9)*ZSIN4+ZFOZQ(10)*ZCOS5+ZFOZQ(11)*ZSIN5)
410  ZOZH(IL)=&
411   &ZFOZH(1)+_TWO_*(ZFOZH(2)*ZCOS1+ZFOZH(3)*ZSIN1+ZFOZH(4)*ZCOS2 &
412   &+ZFOZH(5)*ZSIN2+ZFOZH(6)*ZCOS3+ZFOZH(7)*ZSIN3+ZFOZH(8)&
413   &*ZCOS4+ZFOZH(9)*ZSIN4+ZFOZH(10)*ZCOS5+ZFOZH(11)*ZSIN5)
414  ZOZH(IL)=SQRT(ZOZH(IL))**3
415ENDDO
416
417!     ------------------------------------------------------------------
418
419!       3.     AEROSOLS
420!              --------
421!***
422!       3.1     CALL TO LEGTRI
423
424KCP__RADACA=11
425KDIM_RADACA=66
426!***
427CALL LEGTRI (ZSIN,KCP__RADACA,KDIM_RADACA,ZALP)
428!***
429
430!       3.2     LEGENDRE TRANSFORM FOR AEROSOLS
431!               -------------------------------
432
433DO JMM=1,21
434  ZFAES(JMM) = _ZERO_
435  ZFAEL(JMM) = _ZERO_
436  ZFAEU(JMM) = _ZERO_
437  ZFAED(JMM) = _ZERO_
438ENDDO
439IMM  = 0
440IMNC = 0
441IMNS = 0
442DO JMM=1,11
443  IMM  = IMM+1
444  DO JNN=JMM,11
445    IMNC = IMNC+1
446    ZFAES(IMM) = ZFAES(IMM)+ZALP(IMNC)*RAESC(IMNC)
447    ZFAEL(IMM) = ZFAEL(IMM)+ZALP(IMNC)*RAELC(IMNC)
448    ZFAEU(IMM) = ZFAEU(IMM)+ZALP(IMNC)*RAEUC(IMNC)
449    ZFAED(IMM) = ZFAED(IMM)+ZALP(IMNC)*RAEDC(IMNC)
450  ENDDO
451  IF(JMM /= 1) THEN
452    IMM  = IMM+1
453    DO JNN=JMM,11
454      IMNS = IMNS+1
455      ZFAES(IMM) = ZFAES(IMM)+ZALP(IMNS+11)*RAESS(IMNS)
456      ZFAEL(IMM) = ZFAEL(IMM)+ZALP(IMNS+11)*RAELS(IMNS)
457      ZFAEU(IMM) = ZFAEU(IMM)+ZALP(IMNS+11)*RAEUS(IMNS)
458      ZFAED(IMM) = ZFAED(IMM)+ZALP(IMNS+11)*RAEDS(IMNS)
459    ENDDO
460  ENDIF
461ENDDO
462
463!       3.3     FOURIER TRANSFORM FOR AEROSOLS
464!               ------------------------------
465
466IL = KSHIFT
467DO JL=KIDIA,KFDIA
468  IL = IL+1
469  ZCOS1    = PCLON(JL)
470  ZSIN1    = PSLON(JL)
471  ZCOS2    = ZCOS1*ZCOS1-ZSIN1*ZSIN1
472  ZSIN2    = ZSIN1*ZCOS1+ZCOS1*ZSIN1
473  ZCOS3    = ZCOS2*ZCOS1-ZSIN2*ZSIN1
474  ZSIN3    = ZSIN2*ZCOS1+ZCOS2*ZSIN1
475  ZCOS4    = ZCOS3*ZCOS1-ZSIN3*ZSIN1
476  ZSIN4    = ZSIN3*ZCOS1+ZCOS3*ZSIN1
477  ZCOS5    = ZCOS4*ZCOS1-ZSIN4*ZSIN1
478  ZSIN5    = ZSIN4*ZCOS1+ZCOS4*ZSIN1
479  ZCOS6    = ZCOS5*ZCOS1-ZSIN5*ZSIN1
480  ZSIN6    = ZSIN5*ZCOS1+ZCOS5*ZSIN1
481  ZCOS7    = ZCOS6*ZCOS1-ZSIN6*ZSIN1
482  ZSIN7    = ZSIN6*ZCOS1+ZCOS6*ZSIN1
483  ZCOS8    = ZCOS7*ZCOS1-ZSIN7*ZSIN1
484  ZSIN8    = ZSIN7*ZCOS1+ZCOS7*ZSIN1
485  ZCOS9    = ZCOS8*ZCOS1-ZSIN8*ZSIN1
486  ZSIN9    = ZSIN8*ZCOS1+ZCOS8*ZSIN1
487  ZCOS10   = ZCOS9*ZCOS1-ZSIN9*ZSIN1
488  ZSIN10   = ZSIN9*ZCOS1+ZCOS9*ZSIN1
489  ZAES(IL) = ZFAES(1) + _TWO_*&
490   &( ZFAES(2)*ZCOS1  + ZFAES(3)*ZSIN1  + ZFAES(4)*ZCOS2 &
491   &+ ZFAES(5)*ZSIN2  + ZFAES(6)*ZCOS3  + ZFAES(7)*ZSIN3 &
492   &+ ZFAES(8)*ZCOS4  + ZFAES(9)*ZSIN4  + ZFAES(10)*ZCOS5 &
493   &+ ZFAES(11)*ZSIN5 + ZFAES(12)*ZCOS6 + ZFAES(13)*ZSIN6 &
494   &+ ZFAES(14)*ZCOS7 + ZFAES(15)*ZSIN7 + ZFAES(16)*ZCOS8 &
495   &+ ZFAES(17)*ZSIN8 + ZFAES(18)*ZCOS9 + ZFAES(19)*ZSIN9 &
496   &+ ZFAES(20)*ZCOS10+ ZFAES(21)*ZSIN10                 )
497  ZAEL(IL) = ZFAEL(1) + _TWO_*&
498   &( ZFAEL(2)*ZCOS1  + ZFAEL(3)*ZSIN1  + ZFAEL(4)*ZCOS2 &
499   &+ ZFAEL(5)*ZSIN2  + ZFAEL(6)*ZCOS3  + ZFAEL(7)*ZSIN3 &
500   &+ ZFAEL(8)*ZCOS4  + ZFAEL(9)*ZSIN4  + ZFAEL(10)*ZCOS5 &
501   &+ ZFAEL(11)*ZSIN5 + ZFAEL(12)*ZCOS6 + ZFAEL(13)*ZSIN6 &
502   &+ ZFAEL(14)*ZCOS7 + ZFAEL(15)*ZSIN7 + ZFAEL(16)*ZCOS8 &
503   &+ ZFAEL(17)*ZSIN8 + ZFAEL(18)*ZCOS9 + ZFAEL(19)*ZSIN9 &
504   &+ ZFAEL(20)*ZCOS10+ ZFAEL(21)*ZSIN10                 )
505  ZAEU(IL) = ZFAEU(1) + _TWO_*&
506   &( ZFAEU(2)*ZCOS1  + ZFAEU(3)*ZSIN1  + ZFAEU(4)*ZCOS2 &
507   &+ ZFAEU(5)*ZSIN2  + ZFAEU(6)*ZCOS3  + ZFAEU(7)*ZSIN3 &
508   &+ ZFAEU(8)*ZCOS4  + ZFAEU(9)*ZSIN4  + ZFAEU(10)*ZCOS5 &
509   &+ ZFAEU(11)*ZSIN5 + ZFAEU(12)*ZCOS6 + ZFAEU(13)*ZSIN6 &
510   &+ ZFAEU(14)*ZCOS7 + ZFAEU(15)*ZSIN7 + ZFAEU(16)*ZCOS8 &
511   &+ ZFAEU(17)*ZSIN8 + ZFAEU(18)*ZCOS9 + ZFAEU(19)*ZSIN9 &
512   &+ ZFAEU(20)*ZCOS10+ ZFAEU(21)*ZSIN10                 )
513  ZAED(IL) = ZFAED(1) + _TWO_*&
514   &( ZFAED(2)*ZCOS1  + ZFAED(3)*ZSIN1  + ZFAED(4)*ZCOS2 &
515   &+ ZFAED(5)*ZSIN2  + ZFAED(6)*ZCOS3  + ZFAED(7)*ZSIN3 &
516   &+ ZFAED(8)*ZCOS4  + ZFAED(9)*ZSIN4  + ZFAED(10)*ZCOS5 &
517   &+ ZFAED(11)*ZSIN5 + ZFAED(12)*ZCOS6 + ZFAED(13)*ZSIN6 &
518   &+ ZFAED(14)*ZCOS7 + ZFAED(15)*ZSIN7 + ZFAED(16)*ZCOS8 &
519   &+ ZFAED(17)*ZSIN8 + ZFAED(18)*ZCOS9 + ZFAED(19)*ZSIN9 &
520   &+ ZFAED(20)*ZCOS10+ ZFAED(21)*ZSIN10                 )
521ENDDO
522
523
524!     ------------------------------------------------------------------
525
526!*       4.      VERTICAL DISTRIBUTION
527!*               ---------------------
528
529
530IL=KSHIFT
531DO JL=KIDIA,KFDIA
532  IL=IL+1
533  ZDPO(IL)=PAPRS (JL,1)
534  ZCPHO3=PAPRS (JL,1)**3
535  ZSDPO3=SQRT  (ZCPHO3)
536  IF (LNEWAER) THEN
537    ZAEQSO(IL)= ZAERSS(IL)*CVDAES(1)
538    ZAEQLO(IL)=(ZAEROR(IL)+ZAERSU(IL))*CVDAEL(1)
539    ZAEQUO(IL)= ZAERBC(IL)*CVDAEU(1)
540    ZAEQDO(IL)= ZAERSD(IL)*CVDAED(1)
541  ELSE
542    ZAEQSO(IL)=RCAEOPS*ZAES(IL)*CVDAES(1)
543    ZAEQLO(IL)=RCAEOPL*ZAEL(IL)*CVDAEL(1)
544    ZAEQUO(IL)=RCAEOPU*ZAEU(IL)*CVDAEU(1)
545    ZAEQDO(IL)=RCAEOPD*ZAED(IL)*CVDAED(1)
546  END IF 
547  ZAETRO(IL)=_ONE_
548  ZQOFO(IL)=ZOZQ(IL)*ZSDPO3 / (ZSDPO3 + ZOZH(IL))
549ENDDO
550
551DO JK=1,KLEV
552  IL=KSHIFT
553  IF (KCF == 0) THEN
554    DO JL=KIDIA,KFDIA
555      IL=IL+1
556      ZGRTH(IL)= PTH(JL,JK)/PTH(JL,JK+1)
557    ENDDO
558  ELSEIF (KCF == 1) THEN
559    DO JL=KIDIA,KFDIA
560      IL=IL+1
561      ZGRTH(IL)= PTH(IL,JK)/PTH(IL,JK+1)
562    ENDDO
563  ENDIF
564
565  IL=KSHIFT
566  DO JL=KIDIA,KFDIA
567    IL=IL+1
568    ZDPN(IL)=PAPRS (JL,JK+1)
569    ZCPHN3=PAPRS (JL,JK+1)**3
570    ZSDPN3=SQRT  (ZCPHN3)
571    IF (LNEWAER) THEN
572      ZAEQSN(IL)= ZAERSS(IL)*CVDAES(JK+1)
573      ZAEQLN(IL)=(ZAEROR(IL)+ZAERSU(IL))*CVDAEL(JK+1)
574      ZAEQUN(IL)= ZAERBC(IL)*CVDAEU(JK+1)
575      ZAEQDN(IL)= ZAERSD(IL)*CVDAED(JK+1)
576    ELSE
577      ZAEQSN(IL)=RCAEOPS*ZAES(IL)*CVDAES(JK+1)
578      ZAEQLN(IL)=RCAEOPL*ZAEL(IL)*CVDAEL(JK+1)
579      ZAEQUN(IL)=RCAEOPU*ZAEU(IL)*CVDAEU(JK+1)
580      ZAEQDN(IL)=RCAEOPD*ZAED(IL)*CVDAED(JK+1)
581    END IF 
582
583    IF (_HALF_*(PAPRS(JL,JK)+PAPRS(JL,JK+1)) < 999._JPRB) THEN
584! for models with top above 10hPa
585      ZAETRN(IL)=_ONE_
586      ZAETRO(IL)=_ONE_
587    ELSE
588      ZAETRN(IL)=ZAETRO(IL)*(MIN(_ONE_,     ZGRTH(IL)         ))**RCTRPT
589    ENDIF
590
591    ZAETR=SQRT  (ZAETRN(IL)*ZAETRO(IL))
592    ZQOFN(IL)=ZOZQ(IL)*ZSDPN3/(ZSDPN3+ZOZH(IL))
593    ZDPNMO    =ZDPN(IL)-ZDPO(IL)
594    PAER(IL,1,JK)=(_ONE_-ZAETR)*(RCTRBGA*ZDPNMO+ ZAEQLN(IL)-ZAEQLO(IL))
595    PAER(IL,2,JK)=(_ONE_-ZAETR)*(ZAEQSN(IL)-ZAEQSO(IL))
596    PAER(IL,3,JK)=(_ONE_-ZAETR)*(ZAEQDN(IL)-ZAEQDO(IL))
597    PAER(IL,4,JK)=(_ONE_-ZAETR)*(ZAEQUN(IL)-ZAEQUO(IL))
598!old volc  PAER(IL,JK,5)= ZAETR * RCVOBGA*ZDPNMO
599    PAER(IL,5,JK)=   ZAETR  * ZAERVO(IL) * ZDPNMO
600    PAER(IL,6,JK)=   ZAETR  * RCSTBGA*ZDPNMO
601!old RH dependence         
602!         AADS(IL,JK)=MAX(RCAEADM, (RCAEADK(1)*PAER(IL,1,JK)
603!           + RCAEADK(2)*PAER(IL,2,JK)+RCAEADK(3)*PAER(IL,3,JK))/ZDPNMO)
604    POZON(IL,JK)=ZQOFN(IL)-ZQOFO(IL)
605!**** **************************************************
606!**** **************************************************
607  ENDDO
608  IL=KSHIFT
609  DO JL=KIDIA,KFDIA
610    IL=IL+1
611    ZDPO(IL)=ZDPN(IL)
612    ZQOFO(IL)=ZQOFN(IL)
613
614    ZAEQSO(IL)=ZAEQSN(IL)
615    ZAEQLO(IL)=ZAEQLN(IL)
616    ZAEQUO(IL)=ZAEQUN(IL)
617    ZAEQDO(IL)=ZAEQDN(IL)
618    ZAETRO(IL)=ZAETRN(IL)
619  ENDDO
620 
621 
622!-- diagnostics in case of problem 
623  DO JAER=1,6
624    IL=KSHIFT
625    DO JL=KIDIA,KFDIA
626      IL=IL+1
627      PAER(IL,JAER,JK)=MAX(PAER(IL,JAER,JK),REPAER)
628    END DO
629    itot=il
630  END DO   
631!--
632   
633ENDDO
634
635!     ------------------------------------------------------------------
636
637RETURN
638END SUBROUTINE RADACA
Note: See TracBrowser for help on using the repository browser.