/[lmdze]/trunk/dyn3d/calfis.f
ViewVC logotype

Annotation of /trunk/dyn3d/calfis.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 39 - (hide annotations)
Tue Jan 25 15:11:05 2011 UTC (13 years, 4 months ago) by guez
Original Path: trunk/libf/dyn3d/calfis.f90
File size: 11332 byte(s)
"pi" comes from "nr_util". Removed subroutine "initialize" in module
"comconst".

Copied the content of "fxy_sin.h" into "fxysinus", instead of getting
it from an "include" line. Removed file "fxy_sin.h".

"ps" has rank 2 in "gcm" and "dynetat0".

Assumed-shape for argument "q" of "integrd".

1 guez 3 module calfis_m
2    
3     IMPLICIT NONE
4    
5     contains
6    
7 guez 37 SUBROUTINE calfis(rdayvrai, heure, pucov, pvcov, pteta, q, &
8 guez 10 pmasse, pps, ppk, pphis, pphi, pducov, pdvcov, pdteta, pdq, pw, &
9 guez 37 pdufi, pdvfi, pdhfi, pdqfi, pdpsfi, lafin)
10 guez 3
11 guez 35 ! From dyn3d/calfis.F, version 1.3 2005/05/25 13:10:09
12     ! Authors : P. Le Van, F. Hourdin
13 guez 3
14     ! 1. rearrangement des tableaux et transformation
15     ! variables dynamiques > variables physiques
16     ! 2. calcul des termes physiques
17     ! 3. retransformation des tendances physiques en tendances dynamiques
18    
19     ! remarques:
20     ! ----------
21    
22     ! - les vents sont donnes dans la physique par leurs composantes
23     ! naturelles.
24     ! - la variable thermodynamique de la physique est une variable
25     ! intensive : T
26     ! pour la dynamique on prend T * (preff / p(l)) **kappa
27     ! - les deux seules variables dependant de la geometrie necessaires
28     ! pour la physique sont la latitude pour le rayonnement et
29     ! l'aire de la maille quand on veut integrer une grandeur
30     ! horizontalement.
31    
32     ! Input :
33     ! -------
34     ! pucov covariant zonal velocity
35     ! pvcov covariant meridional velocity
36     ! pteta potential temperature
37     ! pps surface pressure
38     ! pmasse masse d'air dans chaque maille
39     ! pts surface temperature (K)
40     ! callrad clef d'appel au rayonnement
41    
42     ! Output :
43     ! --------
44     ! pdufi tendency for the natural zonal velocity (ms-1)
45     ! pdvfi tendency for the natural meridional velocity
46     ! pdhfi tendency for the potential temperature
47     ! pdtsfi tendency for the surface temperature
48    
49     ! pdtrad radiative tendencies \ both input
50     ! pfluxrad radiative fluxes / and output
51    
52 guez 39 use comconst, only: kappa, cpp, dtphys, g
53 guez 20 use comvert, only: preff
54 guez 3 use comgeom, only: apoln, cu_2d, cv_2d, unsaire_2d, apols, rlonu, rlonv
55 guez 35 use dimens_m, only: iim, jjm, llm, nqmx
56     use dimphy, only: klon
57     use grid_change, only: dyn_phy, gr_fi_dyn
58 guez 18 use iniadvtrac_m, only: niadv
59 guez 39 use nr_util, only: pi
60 guez 3 use physiq_m, only: physiq
61 guez 10 use pressure_var, only: p3d, pls
62 guez 3
63     ! Arguments :
64    
65     LOGICAL, intent(in):: lafin
66     REAL, intent(in):: heure ! heure de la journée en fraction de jour
67    
68 guez 34 REAL pvcov(iim + 1, jjm, llm)
69     REAL pucov(iim + 1, jjm + 1, llm)
70     REAL pteta(iim + 1, jjm + 1, llm)
71     REAL pmasse(iim + 1, jjm + 1, llm)
72 guez 3
73 guez 34 REAL, intent(in):: q(iim + 1, jjm + 1, llm, nqmx)
74 guez 3 ! (mass fractions of advected fields)
75    
76 guez 34 REAL pphis(iim + 1, jjm + 1)
77     REAL pphi(iim + 1, jjm + 1, llm)
78 guez 3
79 guez 34 REAL pdvcov(iim + 1, jjm, llm)
80     REAL pducov(iim + 1, jjm + 1, llm)
81     REAL pdteta(iim + 1, jjm + 1, llm)
82     REAL pdq(iim + 1, jjm + 1, llm, nqmx)
83 guez 3
84 guez 34 REAL pw(iim + 1, jjm + 1, llm)
85 guez 3
86 guez 34 REAL pps(iim + 1, jjm + 1)
87     REAL, intent(in):: ppk(iim + 1, jjm + 1, llm)
88 guez 3
89 guez 34 REAL pdvfi(iim + 1, jjm, llm)
90     REAL pdufi(iim + 1, jjm + 1, llm)
91     REAL pdhfi(iim + 1, jjm + 1, llm)
92     REAL pdqfi(iim + 1, jjm + 1, llm, nqmx)
93     REAL pdpsfi(iim + 1, jjm + 1)
94 guez 3
95     ! Local variables :
96    
97 guez 34 INTEGER i, j, l, ig0, ig, iq, iiq
98 guez 3 REAL zpsrf(klon)
99 guez 34 REAL zplev(klon, llm+1), zplay(klon, llm)
100     REAL zphi(klon, llm), zphis(klon)
101 guez 3
102 guez 35 REAL zufi(klon, llm), v(klon, llm)
103     real zvfi(iim + 1, jjm + 1, llm)
104 guez 34 REAL ztfi(klon, llm) ! temperature
105     real qx(klon, llm, nqmx) ! mass fractions of advected fields
106     REAL pvervel(klon, llm)
107 guez 3
108 guez 34 REAL zdufi(klon, llm), zdvfi(klon, llm)
109     REAL zdtfi(klon, llm), zdqfi(klon, llm, nqmx)
110 guez 3 REAL zdpsrf(klon)
111    
112 guez 35 REAL z1(iim)
113 guez 34 REAL pksurcp(iim + 1, jjm + 1)
114 guez 3
115     ! I. Musat: diagnostic PVteta, Amip2
116     INTEGER, PARAMETER:: ntetaSTD=3
117     REAL:: rtetaSTD(ntetaSTD) = (/350., 380., 405./)
118 guez 34 REAL PVteta(klon, ntetaSTD)
119 guez 3
120 guez 7 REAL, intent(in):: rdayvrai
121 guez 3
122     !-----------------------------------------------------------------------
123    
124     !!print *, "Call sequence information: calfis"
125    
126     ! 1. Initialisations :
127     ! latitude, longitude et aires des mailles pour la physique:
128    
129     ! 40. transformation des variables dynamiques en variables physiques:
130     ! 41. pressions au sol (en Pascals)
131    
132 guez 34 zpsrf(1) = pps(1, 1)
133 guez 3
134     ig0 = 2
135 guez 34 DO j = 2, jjm
136     CALL SCOPY(iim, pps(1, j), 1, zpsrf(ig0), 1)
137 guez 3 ig0 = ig0+iim
138     ENDDO
139    
140 guez 34 zpsrf(klon) = pps(1, jjm + 1)
141 guez 3
142     ! 42. pression intercouches :
143    
144     ! .... zplev definis aux (llm +1) interfaces des couches ....
145     ! .... zplay definis aux (llm) milieux des couches ....
146    
147     ! ... Exner = cp * (p(l) / preff) ** kappa ....
148    
149 guez 10 forall (l = 1: llm+1) zplev(:, l) = pack(p3d(:, :, l), dyn_phy)
150 guez 3
151     ! 43. temperature naturelle (en K) et pressions milieux couches .
152 guez 34 DO l=1, llm
153 guez 10 pksurcp = ppk(:, :, l) / cpp
154     pls(:, :, l) = preff * pksurcp**(1./ kappa)
155     zplay(:, l) = pack(pls(:, :, l), dyn_phy)
156     ztfi(:, l) = pack(pteta(:, :, l) * pksurcp, dyn_phy)
157 guez 3 ENDDO
158    
159     ! 43.bis traceurs
160 guez 34 DO iq=1, nqmx
161 guez 3 iiq=niadv(iq)
162 guez 34 DO l=1, llm
163     qx(1, l, iq) = q(1, 1, l, iiq)
164 guez 3 ig0 = 2
165 guez 34 DO j=2, jjm
166 guez 3 DO i = 1, iim
167 guez 34 qx(ig0, l, iq) = q(i, j, l, iiq)
168 guez 3 ig0 = ig0 + 1
169     ENDDO
170     ENDDO
171 guez 34 qx(ig0, l, iq) = q(1, jjm + 1, l, iiq)
172 guez 3 ENDDO
173     ENDDO
174    
175     ! Geopotentiel calcule par rapport a la surface locale:
176     forall (l = 1:llm) zphi(:, l) = pack(pphi(:, :, l), dyn_phy)
177     zphis = pack(pphis, dyn_phy)
178 guez 34 DO l=1, llm
179     DO ig=1, klon
180     zphi(ig, l)=zphi(ig, l)-zphis(ig)
181 guez 3 ENDDO
182     ENDDO
183    
184 guez 35 ! Calcul de la vitesse verticale (en Pa*m*s ou Kg/s)
185 guez 34 DO l=1, llm
186     pvervel(1, l)=pw(1, 1, l) * g /apoln
187 guez 3 ig0=2
188 guez 34 DO j=2, jjm
189 guez 3 DO i = 1, iim
190 guez 34 pvervel(ig0, l) = pw(i, j, l) * g * unsaire_2d(i, j)
191 guez 3 ig0 = ig0 + 1
192     ENDDO
193     ENDDO
194 guez 34 pvervel(ig0, l)=pw(1, jjm + 1, l) * g /apols
195 guez 3 ENDDO
196    
197     ! 45. champ u:
198    
199 guez 34 DO l=1, llm
200     DO j=2, jjm
201 guez 3 ig0 = 1+(j-2)*iim
202 guez 34 zufi(ig0+1, l)= 0.5 * &
203     (pucov(iim, j, l)/cu_2d(iim, j) + pucov(1, j, l)/cu_2d(1, j))
204     DO i=2, iim
205     zufi(ig0+i, l)= 0.5 * &
206     (pucov(i-1, j, l)/cu_2d(i-1, j) &
207     + pucov(i, j, l)/cu_2d(i, j))
208 guez 3 end DO
209     end DO
210     end DO
211    
212     ! 46.champ v:
213    
214 guez 35 forall (j = 2: jjm, l = 1: llm) zvfi(:iim, j, l)= 0.5 &
215     * (pvcov(:iim, j-1, l) / cv_2d(:iim, j-1) &
216     + pvcov(:iim, j, l) / cv_2d(:iim, j))
217     zvfi(iim + 1, 2:jjm, :) = zvfi(1, 2:jjm, :)
218 guez 3
219 guez 34 ! 47. champs de vents au pôle nord
220 guez 3 ! U = 1 / pi * integrale [ v * cos(long) * d long ]
221     ! V = 1 / pi * integrale [ v * sin(long) * d long ]
222    
223 guez 34 DO l=1, llm
224     z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1, 1, l)/cv_2d(1, 1)
225     DO i=2, iim
226     z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i, 1, l)/cv_2d(i, 1)
227 guez 3 ENDDO
228    
229 guez 35 zufi(1, l) = SUM(COS(rlonv(:iim)) * z1) / pi
230     zvfi(:, 1, l) = SUM(SIN(rlonv(:iim)) * z1) / pi
231 guez 3 ENDDO
232    
233 guez 34 ! 48. champs de vents au pôle sud:
234 guez 3 ! U = 1 / pi * integrale [ v * cos(long) * d long ]
235     ! V = 1 / pi * integrale [ v * sin(long) * d long ]
236    
237 guez 34 DO l=1, llm
238     z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1, jjm, l) &
239     /cv_2d(1, jjm)
240     DO i=2, iim
241     z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i, jjm, l)/cv_2d(i, jjm)
242 guez 3 ENDDO
243    
244 guez 35 zufi(klon, l) = SUM(COS(rlonv(:iim)) * z1) / pi
245     zvfi(:, jjm + 1, l) = SUM(SIN(rlonv(:iim)) * z1) / pi
246     ENDDO
247 guez 3
248 guez 35 forall(l= 1: llm) v(:, l) = pack(zvfi(:, :, l), dyn_phy)
249 guez 3
250     !IM calcul PV a teta=350, 380, 405K
251 guez 35 CALL PVtheta(klon, llm, pucov, pvcov, pteta, ztfi, zplay, zplev, &
252 guez 34 ntetaSTD, rtetaSTD, PVteta)
253 guez 3
254 guez 35 ! Appel de la physique :
255     CALL physiq(lafin, rdayvrai, heure, dtphys, zplev, zplay, zphi, &
256     zphis, zufi, v, ztfi, qx, pvervel, zdufi, zdvfi, &
257     zdtfi, zdqfi, zdpsrf, pducov, PVteta) ! diagnostic PVteta, Amip2
258 guez 3
259     ! transformation des tendances physiques en tendances dynamiques:
260    
261     ! tendance sur la pression :
262    
263     pdpsfi = gr_fi_dyn(zdpsrf)
264    
265     ! 62. enthalpie potentielle
266    
267 guez 34 DO l=1, llm
268 guez 3
269 guez 34 DO i=1, iim + 1
270     pdhfi(i, 1, l) = cpp * zdtfi(1, l) / ppk(i, 1 , l)
271     pdhfi(i, jjm + 1, l) = cpp * zdtfi(klon, l)/ ppk(i, jjm + 1, l)
272 guez 3 ENDDO
273    
274 guez 34 DO j=2, jjm
275 guez 3 ig0=1+(j-2)*iim
276 guez 34 DO i=1, iim
277     pdhfi(i, j, l) = cpp * zdtfi(ig0+i, l) / ppk(i, j, l)
278 guez 3 ENDDO
279 guez 34 pdhfi(iim + 1, j, l) = pdhfi(1, j, l)
280 guez 3 ENDDO
281    
282     ENDDO
283    
284     ! 62. humidite specifique
285    
286 guez 34 DO iq=1, nqmx
287     DO l=1, llm
288     DO i=1, iim + 1
289     pdqfi(i, 1, l, iq) = zdqfi(1, l, iq)
290     pdqfi(i, jjm + 1, l, iq) = zdqfi(klon, l, iq)
291 guez 3 ENDDO
292 guez 34 DO j=2, jjm
293 guez 3 ig0=1+(j-2)*iim
294 guez 34 DO i=1, iim
295     pdqfi(i, j, l, iq) = zdqfi(ig0+i, l, iq)
296 guez 3 ENDDO
297 guez 34 pdqfi(iim + 1, j, l, iq) = pdqfi(1, j, l, iq)
298 guez 3 ENDDO
299     ENDDO
300     ENDDO
301    
302     ! 63. traceurs
303    
304     ! initialisation des tendances
305     pdqfi=0.
306    
307 guez 34 DO iq=1, nqmx
308 guez 3 iiq=niadv(iq)
309 guez 34 DO l=1, llm
310     DO i=1, iim + 1
311     pdqfi(i, 1, l, iiq) = zdqfi(1, l, iq)
312     pdqfi(i, jjm + 1, l, iiq) = zdqfi(klon, l, iq)
313 guez 3 ENDDO
314 guez 34 DO j=2, jjm
315 guez 3 ig0=1+(j-2)*iim
316 guez 34 DO i=1, iim
317     pdqfi(i, j, l, iiq) = zdqfi(ig0+i, l, iq)
318 guez 3 ENDDO
319 guez 34 pdqfi(iim + 1, j, l, iiq) = pdqfi(1, j, l, iq)
320 guez 3 ENDDO
321     ENDDO
322     ENDDO
323    
324     ! 65. champ u:
325    
326 guez 34 DO l=1, llm
327 guez 3
328 guez 34 DO i=1, iim + 1
329     pdufi(i, 1, l) = 0.
330     pdufi(i, jjm + 1, l) = 0.
331 guez 3 ENDDO
332    
333 guez 34 DO j=2, jjm
334 guez 3 ig0=1+(j-2)*iim
335 guez 34 DO i=1, iim-1
336     pdufi(i, j, l)= &
337     0.5*(zdufi(ig0+i, l)+zdufi(ig0+i+1, l))*cu_2d(i, j)
338 guez 3 ENDDO
339 guez 34 pdufi(iim, j, l)= &
340     0.5*(zdufi(ig0+1, l)+zdufi(ig0+iim, l))*cu_2d(iim, j)
341     pdufi(iim + 1, j, l)=pdufi(1, j, l)
342 guez 3 ENDDO
343    
344     ENDDO
345    
346     ! 67. champ v:
347    
348 guez 34 DO l=1, llm
349 guez 3
350 guez 34 DO j=2, jjm-1
351 guez 3 ig0=1+(j-2)*iim
352 guez 34 DO i=1, iim
353     pdvfi(i, j, l)= &
354     0.5*(zdvfi(ig0+i, l)+zdvfi(ig0+i+iim, l))*cv_2d(i, j)
355 guez 3 ENDDO
356 guez 34 pdvfi(iim + 1, j, l) = pdvfi(1, j, l)
357 guez 3 ENDDO
358     ENDDO
359    
360     ! 68. champ v pres des poles:
361     ! v = U * cos(long) + V * SIN(long)
362    
363 guez 34 DO l=1, llm
364     DO i=1, iim
365     pdvfi(i, 1, l)= &
366     zdufi(1, l)*COS(rlonv(i))+zdvfi(1, l)*SIN(rlonv(i))
367     pdvfi(i, jjm, l)=zdufi(klon, l)*COS(rlonv(i)) &
368     +zdvfi(klon, l)*SIN(rlonv(i))
369     pdvfi(i, 1, l)= &
370     0.5*(pdvfi(i, 1, l)+zdvfi(i+1, l))*cv_2d(i, 1)
371     pdvfi(i, jjm, l)= &
372     0.5*(pdvfi(i, jjm, l)+zdvfi(klon-iim-1+i, l))*cv_2d(i, jjm)
373 guez 3 ENDDO
374    
375 guez 34 pdvfi(iim + 1, 1, l) = pdvfi(1, 1, l)
376     pdvfi(iim + 1, jjm, l)= pdvfi(1, jjm, l)
377 guez 3 ENDDO
378    
379     END SUBROUTINE calfis
380    
381     end module calfis_m

  ViewVC Help
Powered by ViewVC 1.1.21