source: trunk/LMDZ.MARS/libf/phymars/flusv.F @ 1266

Last change on this file since 1266 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

  • Optimized computations in paramfoto_compact (twice less dlog10 calculations)
  • Checked consistency before/after modification in debug mode
  • Checked performance is not impacted (same as before)
File size: 13.7 KB
Line 
1      SUBROUTINE flusv(KDLON,nsf,n,omega,g,tau,emis,bh,bsol,fah,fdh)
2      use dimradmars_mod, only: ndlo2, ndlon, nflev
3      IMPLICIT NONE
4c.......................................................................
5c
6c  calcul des flux ascendant et descendant aux interfaces entre n couches
7c  * dans l'infrarouge
8c  * B est une fonction lineaire de $\tau$ a l'interieur de chaque couche
9c  * le B du sol peut etre different de celui qui correspond au profil
10c    de la n-ieme couche
11c  * l'hypothese est une hypothese a deux flux isotropes sur chaque
12c    hemisphere ("hemispheric constant") + "source function technique"
13c    (voir Toon et al. 1988)
14c  * le flux descendant en haut de l'atmosphere est nul
15c  * les couches sont numerotees du haut de l'atmosphere vers le sol
16c
17c  in :   * KDLON      ---> dimension de vectorisation
18c         * nsf        ---> nsf=0 ==> "hemispheric constant"
19c                           nsf>0 ==> "hemispheric constant" + "source function"
20c         * n          ---> nombre de couches
21c         * omega(i)   ---> single scattering albedo pour la i-eme couche
22c         * g(i)       ---> asymmetry parameter pour la i-eme couche
23c         * tau(i)     ---> epaisseur optique de la i-eme couche
24c         * emis       ---> emissivite du sol
25c         * bh(i)      ---> luminance du corps noir en haut de la i-eme
26c                           couche, bh(n+1) pour la valeur au sol qui
27c                           correspond au profil de la n-ieme couche
28c         * bsol       ---> luminance du corps noir au sol
29c
30c  out :  * fah(i)     ---> flux ascendant en haut de la i-eme couche,
31c                           fah(n+1) pour le sol
32c         * fdh(i)     ---> flux descendant en haut de la i-eme couche,
33c                           fdh(n+1) pour le sol
34c
35c.......................................................................
36c  include des dimensions locales
37c
38#include "dimensions.h"
39c.......................................................................
40c  declaration des arguments
41c
42      INTEGER KDLON,nsf,n
43      REAL omega(NDLO2,n),g(NDLO2,n),tau(NDLO2,n),emis(NDLO2)
44     &,bh(NDLO2,n+1),bsol(NDLO2),fah(NDLO2,n+1),fdh(NDLO2,n+1)
45c.......................................................................
46c  declaration des variables locales
47c
48      REAL pi
49      PARAMETER (pi=3.141592653589793E+0)
50      INTEGER iv,i,j
51      REAL beta,gama1,gama2,amu1,grgama,b0,b1
52      REAL a(NDLON,4*nflev),b(NDLON,4*nflev)
53     &    ,d(NDLON,4*nflev),e(NDLON,4*nflev)
54     &    ,y(NDLON,4*nflev)
55     &    ,alambda(NDLON,2*nflev)
56     &    ,e1(NDLON,2*nflev),e2(NDLON,2*nflev)
57     &    ,e3(NDLON,2*nflev),e4(NDLON,2*nflev)
58     &    ,cah(NDLON,2*nflev),cab(NDLON,2*nflev)
59     &    ,cdh(NDLON,2*nflev),cdb(NDLON,2*nflev)
60      REAL grg(NDLON,2*nflev),grh(NDLON,2*nflev)
61     &    ,grj(NDLON,2*nflev),grk(NDLON,2*nflev)
62     &    ,alpha1(NDLON,2*nflev),alpha2(NDLON,2*nflev)
63     &    ,sigma1(NDLON,2*nflev),sigma2(NDLON,2*nflev)
64      INTEGER nq
65      PARAMETER (nq=8)
66      REAL x(nq),w(nq),gri(NDLON,nq)
67      DATA x/1.9855071751231860E-2 , 0.1016667612931866E+0
68     &     , 0.2372337950418355E+0 , 0.4082826787521751E+0
69     &     , 0.5917173212478250E+0 , 0.7627662049581645E+0
70     &     , 0.8983332387068134E+0 , 0.9801449282487682E+0/
71      DATA w/5.0614268145185310E-2 , 0.1111905172266872E+0
72     &     , 0.1568533229389437E+0 , 0.1813418916891810E+0
73     &     , 0.1813418916891810E+0 , 0.1568533229389437E+0
74     &     , 0.1111905172266872E+0 , 5.0614268145185310E-2/
75c.......................................................................
76c
77c.......................................................................
78      do 10001 i=1,n
79                                                   do 99999 iv=1,KDLON
80      beta=(1.E+0-g(iv,i))/2.E+0
81      gama1=2.E+0*(1.E+0-omega(iv,i)*(1.E+0-beta))
82      gama2=2.E+0*omega(iv,i)*beta
83      amu1=5.E-1
84      alambda(iv,i)=sqrt(gama1**2-gama2**2)
85      grgama=(gama1-alambda(iv,i))/gama2
86c
87c ici une petite bidouille : si l'epaisseur optique d'une couche
88c est trop faible, $dB \over d\tau$ devient tres grand et le schema
89c s'ecroule. dans ces cas, on fait l'hypothese de couche isotherme.
90c
91      if (tau(iv,i).gt.1.E-3) then
92      b0=bh(iv,i)
93      b1=(bh(iv,i+1)-b0)/tau(iv,i)
94      else
95      b0=(bh(iv,i)+bh(iv,i+1))/2.E+0
96      b1=0.E+0
97      endif
98c
99      e1(iv,i)=1.E+0+grgama*exp(-alambda(iv,i)*tau(iv,i))
100      e2(iv,i)=1.E+0-grgama*exp(-alambda(iv,i)*tau(iv,i))
101      e3(iv,i)=grgama+exp(-alambda(iv,i)*tau(iv,i))
102      e4(iv,i)=grgama-exp(-alambda(iv,i)*tau(iv,i))
103      cah(iv,i)=2.E+0*pi*amu1*(b0+b1/(gama1+gama2))
104      cab(iv,i)=2.E+0*pi*amu1*(b0+b1*(tau(iv,i)+1.E+0/(gama1+gama2)))
105      cdh(iv,i)=2.E+0*pi*amu1*(b0-b1/(gama1+gama2))
106      cdb(iv,i)=2.E+0*pi*amu1*(b0+b1*(tau(iv,i)-1.E+0/(gama1+gama2))) 
107c
108      grg(iv,i)=(1.E+0/amu1-alambda(iv,i))
109      grh(iv,i)=grgama*(alambda(iv,i)+1.E+0/amu1)
110      grj(iv,i)=grh(iv,i)
111      grk(iv,i)=grg(iv,i)
112      alpha1(iv,i)=2.E+0*pi*(b0+b1*(1.E+0/(gama1+gama2)-amu1))
113      alpha2(iv,i)=2.E+0*pi*b1
114      sigma1(iv,i)=2.E+0*pi*(b0-b1*(1.E+0/(gama1+gama2)-amu1))
115      sigma2(iv,i)=alpha2(iv,i)
116c
11799999                                              continue
11810001 continue
119c.......................................................................
120                                                   do 99998 iv=1,KDLON
121      a(iv,1)=0.E+0
122      b(iv,1)=e1(iv,1)
123      d(iv,1)=-e2(iv,1)
124      e(iv,1)=-cdh(iv,1)
12599998                                              continue
126c
127      do 10002 i=1,n-1
128      j=2*i+1
129                                                   do 99997 iv=1,KDLON
130      a(iv,j)=e2(iv,i)*e3(iv,i)-e4(iv,i)*e1(iv,i)
131      b(iv,j)=e1(iv,i)*e1(iv,i+1)-e3(iv,i)*e3(iv,i+1)
132      d(iv,j)=e3(iv,i)*e4(iv,i+1)-e1(iv,i)*e2(iv,i+1)
133      e(iv,j)=e3(iv,i)*(cah(iv,i+1)-cab(iv,i))
134     &       +e1(iv,i)*(cdb(iv,i)-cdh(iv,i+1))
13599997                                              continue
13610002 continue
137c
138      do 10003 i=1,n-1
139      j=2*i
140                                                   do 99996 iv=1,KDLON
141      a(iv,j)=e2(iv,i+1)*e1(iv,i)-e3(iv,i)*e4(iv,i+1)
142      b(iv,j)=e2(iv,i)*e2(iv,i+1)-e4(iv,i)*e4(iv,i+1)
143      d(iv,j)=e1(iv,i+1)*e4(iv,i+1)-e2(iv,i+1)*e3(iv,i+1)
144      e(iv,j)=e2(iv,i+1)*(cah(iv,i+1)-cab(iv,i))
145     &       +e4(iv,i+1)*(cdb(iv,i)-cdh(iv,i+1))
14699996                                              continue
14710003 continue
148c     
149      j=2*n
150                                                   do 99995 iv=1,KDLON
151      a(iv,j)=e1(iv,n)-(1.E+0-emis(iv))*e3(iv,n)
152      b(iv,j)=e2(iv,n)-(1.E+0-emis(iv))*e4(iv,n)
153      d(iv,j)=0.E+0
154      e(iv,j)=emis(iv)*pi*bsol(iv)-cab(iv,n)+(1.E+0-emis(iv))*cdb(iv,n)
15599995                                              continue
156c.......................................................................
157      call sys3v(KDLON,2*n,a,b,d,e,y)
158c.......................................................................
159      do 10004 i=1,n
160                                                   do 99994 iv=1,KDLON
161      grg(iv,i)=grg(iv,i)*(y(iv,2*i-1)+y(iv,2*i))
162      grh(iv,i)=grh(iv,i)*(y(iv,2*i-1)-y(iv,2*i))
163      grj(iv,i)=grj(iv,i)*(y(iv,2*i-1)+y(iv,2*i))
164      grk(iv,i)=grk(iv,i)*(y(iv,2*i-1)-y(iv,2*i))
16599994                                              continue
16610004 continue
167c.......................................................................
168c les valeurs de flux "hemispheric constant"
169c
170      IF (nsf.eq.0) THEN
171      do 10005 i=1,n
172                                                   do 99993 iv=1,KDLON
173      fah(iv,i)=e3(iv,i)*y(iv,2*i-1)-e4(iv,i)*y(iv,2*i)+cah(iv,i)
174      fdh(iv,i)=e1(iv,i)*y(iv,2*i-1)-e2(iv,i)*y(iv,2*i)+cdh(iv,i)
17599993                                              continue
17610005 continue
177                                                   do 99992 iv=1,KDLON
178      fah(iv,n+1)=e1(iv,n)*y(iv,2*n-1)+e2(iv,n)*y(iv,2*n)+cab(iv,n)
179      fdh(iv,n+1)=e3(iv,n)*y(iv,2*n-1)+e4(iv,n)*y(iv,2*n)+cdb(iv,n)
18099992                                              continue
181      GOTO 10100
182      ENDIF
183c.......................................................................
184c passage a "source function"
185c
186c on applique une quadrature de dimension nq
187c x est le vecteur des \mu de la quadrature
188c w est le vecteur des poids correspondants
189c nq est fixe en parameter
190c x et w sont fixes en data
191c
192c.......................................................................
193c on part d'en haut et on descent selon les nq angles pour calculer
194c tous les flux descendants
195c
196      do 10006 j=1,nq
197                                                   do 99991 iv=1,KDLON
198      gri(iv,j)=0.E+0
19999991                                              continue
20010006 continue
201                                                   do 99990 iv=1,KDLON
202      fdh(iv,1)=0.E+0
20399990                                              continue
204      do 10007 i=1,n
205      do 10008 j=1,nq
206                                                   do 99989 iv=1,KDLON
207      gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j))
208     &+grj(iv,i)/(alambda(iv,i)*x(j)+1.E+0)
209     &*(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j))))
210     &+grk(iv,i)/(alambda(iv,i)*x(j)-1.E+0)
211     &*(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i)))
212     &+sigma1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j)))
213     &+sigma2(iv,i)*(x(j)*exp(-tau(iv,i)/x(j))+tau(iv,i)-x(j))
21499989                                              continue
21510008 continue
216                                                   do 99988 iv=1,KDLON
217      fdh(iv,i+1)=0.E+0
21899988                                              continue
219      do 10009 j=1,nq
220                                                   do 99987 iv=1,KDLON
221      fdh(iv,i+1)=fdh(iv,i+1)+w(j)*x(j)*gri(iv,j)
22299987                                              continue
22310009 continue
22410007 continue
225c.......................................................................
226c on applique la condition de reflexion a sol
227c
228                                                   do 99986 iv=1,KDLON
229      fah(iv,n+1)=(1.E+0-emis(iv))*fdh(iv,n+1)+pi*emis(iv)*bsol(iv)
23099986                                              continue
231      do 10010 j=1,nq
232                                                   do 99985 iv=1,KDLON
233      gri(iv,j)=2.E+0*fah(iv,n+1)
23499985                                              continue
23510010 continue
236c.......................................................................
237c on remonte pour calculer tous les flux ascendants
238c 
239      do 10011 i=n,1,-1
240      do 10012 j=1,nq
241                                                   do 99984 iv=1,KDLON
242      gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j))
243     &+grg(iv,i)/(alambda(iv,i)*x(j)-1.E+0)
244     &*(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i)))
245     &+grh(iv,i)/(alambda(iv,i)*x(j)+1.E+0)
246     &*(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j))))
247     &+alpha1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j)))
248     &+alpha2(iv,i)*(x(j)-(tau(iv,i)+x(j))*exp(-tau(iv,i)/x(j)))
24999984                                              continue
25010012 continue
251                                                   do 99983 iv=1,KDLON
252      fah(iv,i)=0.E+0
25399983                                              continue
254      do 10013 j=1,nq
255                                                   do 99982 iv=1,KDLON
256      fah(iv,i)=fah(iv,i)+w(j)*x(j)*gri(iv,j)
25799982                                              continue
25810013 continue
25910011 continue
260c.......................................................................
26110100 continue     
262c.......................................................................
263c
264c.......................................................................
265      RETURN
266      END
267
268c ***************************************************************
269
270
271      SUBROUTINE sys3v(KDLON,n,a,b,d,e,y)
272      use dimradmars_mod, only: ndlon, ndlo2, nflev
273      IMPLICIT NONE
274c.......................................................................
275c
276c  inversion d'un systeme lineaire tridiagonal
277c
278c  |  b1  d1               |   | y1 |   | e1 |
279c  |  a2  b2  d2           |   | y2 |   | e2 |
280c  |      a3  b3  d3       | * | y3 | = | e3 |
281c  |             ....      |   |    |   |    |
282c  |              an  bn   |   | yn |   | en |
283c
284c  in :   * KDLON        --> dimension de vectorisation
285c         * n            --> dimension du systeme
286c         * a,b,d,e      --> voir dessin
287c
288c  out :  * y            --> voir dessin
289c
290c.......................................................................
291c  declaration des arguments
292c
293      INTEGER KDLON,n
294      REAL a(NDLO2,n),b(NDLO2,n),d(NDLO2,n),e(NDLO2,n),y(NDLO2,n)
295c.......................................................................
296c  declaration des variables locales
297c
298      INTEGER iv,i
299      REAL as(NDLON,4*nflev),ds(NDLON,4*nflev)
300     &    ,x(NDLON,4*nflev)
301c.......................................................................
302c
303c.......................................................................
304                                                   do 99999 iv=1,KDLON
305      as(iv,n)=a(iv,n)/b(iv,n)
306      ds(iv,n)=e(iv,n)/b(iv,n)
30799999                                              continue
308      do 10001 i=n-1,1,-1
309                                                   do 99998 iv=1,KDLON
310      x(iv,i)=1.E+0/(b(iv,i)-d(iv,i)*as(iv,i+1))
311      as(iv,i)=a(iv,i)*x(iv,i)
312      ds(iv,i)=(e(iv,i)-d(iv,i)*ds(iv,i+1))*x(iv,i)
31399998                                              continue
31410001 continue
315                                                   do 99997 iv=1,KDLON
316      y(iv,1)=ds(iv,1)
31799997                                              continue
318      do 10002 i=2,n
319                                                   do 99996 iv=1,KDLON
320      y(iv,i)=ds(iv,i)-as(iv,i)*y(iv,i-1)
32199996                                              continue
32210002 continue
323c.......................................................................
324c
325c.......................................................................
326      RETURN
327      END
Note: See TracBrowser for help on using the repository browser.