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

Contents of /trunk/dyn3d/calfis.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 79 - (show annotations)
Fri Feb 28 17:52:47 2014 UTC (10 years, 3 months ago) by guez
Original Path: trunk/dyn3d/calfis.f90
File size: 9981 byte(s)
Moved procedure iniconst inside module comconst. Removed useless
variables of module comconst: im, jm, lllm, imp1, jmp1, lllmm1,
lllmp1, lcl, cotot, unsim. Move definition of dtvr that was in
dynetat0 and etat0 to iniconst. Moved comparison of dtvr from day_step
and start.nc that was in gcm to dynetat0. Moved call to disvert out of
iniconst. Moved call to iniconst in gcm before call to dynetat0.

Removed unused argument pvteta of physiq (not used either in LMDZ).

1 module calfis_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE calfis(rdayvrai, time, ucov, vcov, teta, q, ps, pk, phis, phi, &
8 dudyn, dv, w, dufi, dvfi, dtetafi, dqfi, dpfi, lafin)
9
10 ! From dyn3d/calfis.F, version 1.3 2005/05/25 13:10:09
11 ! Authors: P. Le Van, F. Hourdin
12
13 ! 1. Réarrangement des tableaux et transformation des variables
14 ! dynamiques en variables physiques
15
16 ! 2. Calcul des termes physiques
17 ! 3. Retransformation des tendances physiques en tendances dynamiques
18
19 ! Remarques:
20
21 ! - Les vents sont donnés dans la physique par leurs composantes
22 ! naturelles.
23
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
28 ! - Les deux seules variables dépendant de la géométrie
29 ! nécessaires pour la physique sont la latitude pour le
30 ! rayonnement et l'aire de la maille quand on veut intégrer une
31 ! grandeur horizontalement.
32
33 use comconst, only: kappa, cpp, dtphys, g
34 use comgeom, only: apoln, cu_2d, cv_2d, unsaire_2d, apols, rlonu, rlonv
35 use dimens_m, only: iim, jjm, llm, nqmx
36 use dimphy, only: klon
37 use disvert_m, only: preff
38 use grid_change, only: dyn_phy, gr_fi_dyn
39 use iniadvtrac_m, only: niadv
40 use nr_util, only: pi
41 use physiq_m, only: physiq
42 use pressure_var, only: p3d, pls
43
44 ! Arguments :
45
46 ! Output :
47 ! dvfi tendency for the natural meridional velocity
48 ! dtetafi tendency for the potential temperature
49 ! pdtsfi tendency for the surface temperature
50
51 ! pdtrad radiative tendencies \ input and output
52 ! pfluxrad radiative fluxes / input and output
53
54 REAL, intent(in):: rdayvrai
55 REAL, intent(in):: time ! heure de la journée en fraction de jour
56 REAL, intent(in):: ucov(iim + 1, jjm + 1, llm)
57 ! ucov covariant zonal velocity
58 REAL, intent(in):: vcov(iim + 1, jjm, llm)
59 ! vcov covariant meridional velocity
60 REAL, intent(in):: teta(iim + 1, jjm + 1, llm)
61 ! teta potential temperature
62
63 REAL, intent(in):: q(iim + 1, jjm + 1, llm, nqmx)
64 ! (mass fractions of advected fields)
65
66 REAL, intent(in):: ps(iim + 1, jjm + 1)
67 ! ps surface pressure
68 REAL, intent(in):: pk(iim + 1, jjm + 1, llm)
69 REAL, intent(in):: phis(iim + 1, jjm + 1)
70 REAL, intent(in):: phi(iim + 1, jjm + 1, llm)
71 REAL dudyn(iim + 1, jjm + 1, llm)
72 REAL dv(iim + 1, jjm, llm)
73 REAL, intent(in):: w(iim + 1, jjm + 1, llm)
74
75 REAL, intent(out):: dufi(iim + 1, jjm + 1, llm)
76 ! tendency for the covariant zonal velocity (m2 s-2)
77
78 REAL dvfi(iim + 1, jjm, llm)
79 REAL, intent(out):: dtetafi(iim + 1, jjm + 1, llm)
80 REAL dqfi(iim + 1, jjm + 1, llm, nqmx)
81 REAL dpfi(iim + 1, jjm + 1)
82 LOGICAL, intent(in):: lafin
83
84 ! Local variables :
85
86 INTEGER i, j, l, ig0, ig, iq, iiq
87 REAL zpsrf(klon)
88 REAL paprs(klon, llm+1), play(klon, llm)
89 REAL pphi(klon, llm), pphis(klon)
90
91 REAL u(klon, llm), v(klon, llm)
92 real zvfi(iim + 1, jjm + 1, llm)
93 REAL t(klon, llm) ! temperature
94 real qx(klon, llm, nqmx) ! mass fractions of advected fields
95 REAL omega(klon, llm)
96
97 REAL d_u(klon, llm), d_v(klon, llm) ! tendances physiques du vent (m s-2)
98 REAL d_t(klon, llm), d_qx(klon, llm, nqmx)
99 REAL d_ps(klon)
100
101 REAL z1(iim)
102 REAL pksurcp(iim + 1, jjm + 1)
103
104 ! I. Musat: diagnostic PVteta, Amip2
105 INTEGER, PARAMETER:: ntetaSTD=3
106 REAL:: rtetaSTD(ntetaSTD) = (/350., 380., 405./)
107 REAL PVteta(klon, ntetaSTD)
108
109 !-----------------------------------------------------------------------
110
111 !!print *, "Call sequence information: calfis"
112
113 ! 1. Initialisations :
114 ! latitude, longitude et aires des mailles pour la physique:
115
116 ! 40. transformation des variables dynamiques en variables physiques:
117 ! 41. pressions au sol (en Pascals)
118
119 zpsrf(1) = ps(1, 1)
120
121 ig0 = 2
122 DO j = 2, jjm
123 CALL SCOPY(iim, ps(1, j), 1, zpsrf(ig0), 1)
124 ig0 = ig0+iim
125 ENDDO
126
127 zpsrf(klon) = ps(1, jjm + 1)
128
129 ! 42. pression intercouches :
130
131 ! paprs defini aux (llm +1) interfaces des couches
132 ! play defini aux (llm) milieux des couches
133
134 ! Exner = cp * (p(l) / preff) ** kappa
135
136 forall (l = 1: llm+1) paprs(:, l) = pack(p3d(:, :, l), dyn_phy)
137
138 ! 43. temperature naturelle (en K) et pressions milieux couches
139 DO l=1, llm
140 pksurcp = pk(:, :, l) / cpp
141 pls(:, :, l) = preff * pksurcp**(1./ kappa)
142 play(:, l) = pack(pls(:, :, l), dyn_phy)
143 t(:, l) = pack(teta(:, :, l) * pksurcp, dyn_phy)
144 ENDDO
145
146 ! 43.bis traceurs
147 DO iq=1, nqmx
148 iiq=niadv(iq)
149 DO l=1, llm
150 qx(1, l, iq) = q(1, 1, l, iiq)
151 ig0 = 2
152 DO j=2, jjm
153 DO i = 1, iim
154 qx(ig0, l, iq) = q(i, j, l, iiq)
155 ig0 = ig0 + 1
156 ENDDO
157 ENDDO
158 qx(ig0, l, iq) = q(1, jjm + 1, l, iiq)
159 ENDDO
160 ENDDO
161
162 ! Geopotentiel calcule par rapport a la surface locale:
163 forall (l = 1:llm) pphi(:, l) = pack(phi(:, :, l), dyn_phy)
164 pphis = pack(phis, dyn_phy)
165 forall (l = 1:llm) pphi(:, l)=pphi(:, l) - pphis
166
167 ! Calcul de la vitesse verticale (en Pa*m*s ou Kg/s)
168 DO l=1, llm
169 omega(1, l)=w(1, 1, l) * g /apoln
170 ig0=2
171 DO j=2, jjm
172 DO i = 1, iim
173 omega(ig0, l) = w(i, j, l) * g * unsaire_2d(i, j)
174 ig0 = ig0 + 1
175 ENDDO
176 ENDDO
177 omega(ig0, l)=w(1, jjm + 1, l) * g /apols
178 ENDDO
179
180 ! 45. champ u:
181
182 DO l=1, llm
183 DO j=2, jjm
184 ig0 = 1+(j-2)*iim
185 u(ig0+1, l)= 0.5 &
186 * (ucov(iim, j, l) / cu_2d(iim, j) + ucov(1, j, l) / cu_2d(1, j))
187 DO i=2, iim
188 u(ig0+i, l)= 0.5 * (ucov(i-1, j, l)/cu_2d(i-1, j) &
189 + ucov(i, j, l)/cu_2d(i, j))
190 end DO
191 end DO
192 end DO
193
194 ! 46.champ v:
195
196 forall (j = 2: jjm, l = 1: llm) zvfi(:iim, j, l)= 0.5 &
197 * (vcov(:iim, j-1, l) / cv_2d(:iim, j-1) &
198 + vcov(:iim, j, l) / cv_2d(:iim, j))
199 zvfi(iim + 1, 2:jjm, :) = zvfi(1, 2:jjm, :)
200
201 ! 47. champs de vents au pôle nord
202 ! U = 1 / pi * integrale [ v * cos(long) * d long ]
203 ! V = 1 / pi * integrale [ v * sin(long) * d long ]
204
205 DO l=1, llm
206 z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*vcov(1, 1, l)/cv_2d(1, 1)
207 DO i=2, iim
208 z1(i) =(rlonu(i)-rlonu(i-1))*vcov(i, 1, l)/cv_2d(i, 1)
209 ENDDO
210
211 u(1, l) = SUM(COS(rlonv(:iim)) * z1) / pi
212 zvfi(:, 1, l) = SUM(SIN(rlonv(:iim)) * z1) / pi
213 ENDDO
214
215 ! 48. champs de vents au pôle sud:
216 ! U = 1 / pi * integrale [ v * cos(long) * d long ]
217 ! V = 1 / pi * integrale [ v * sin(long) * d long ]
218
219 DO l=1, llm
220 z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*vcov(1, jjm, l) &
221 /cv_2d(1, jjm)
222 DO i=2, iim
223 z1(i) =(rlonu(i)-rlonu(i-1))*vcov(i, jjm, l)/cv_2d(i, jjm)
224 ENDDO
225
226 u(klon, l) = SUM(COS(rlonv(:iim)) * z1) / pi
227 zvfi(:, jjm + 1, l) = SUM(SIN(rlonv(:iim)) * z1) / pi
228 ENDDO
229
230 forall(l= 1: llm) v(:, l) = pack(zvfi(:, :, l), dyn_phy)
231
232 ! Compute potential vorticity at theta = 350, 380 and 405 K:
233 CALL PVtheta(klon, llm, ucov, vcov, teta, t, play, paprs, ntetaSTD, &
234 rtetaSTD, PVteta)
235
236 ! Appel de la physique :
237 CALL physiq(lafin, rdayvrai, time, dtphys, paprs, play, pphi, pphis, u, &
238 v, t, qx, omega, d_u, d_v, d_t, d_qx, d_ps, dudyn)
239
240 ! transformation des tendances physiques en tendances dynamiques:
241
242 ! tendance sur la pression :
243
244 dpfi = gr_fi_dyn(d_ps)
245
246 ! 62. enthalpie potentielle
247 do l=1, llm
248 dtetafi(:, :, l) = cpp * gr_fi_dyn(d_t(:, l)) / pk(:, :, l)
249 end do
250
251 ! 62. humidite specifique
252
253 DO iq=1, nqmx
254 DO l=1, llm
255 DO i=1, iim + 1
256 dqfi(i, 1, l, iq) = d_qx(1, l, iq)
257 dqfi(i, jjm + 1, l, iq) = d_qx(klon, l, iq)
258 ENDDO
259 DO j=2, jjm
260 ig0=1+(j-2)*iim
261 DO i=1, iim
262 dqfi(i, j, l, iq) = d_qx(ig0+i, l, iq)
263 ENDDO
264 dqfi(iim + 1, j, l, iq) = dqfi(1, j, l, iq)
265 ENDDO
266 ENDDO
267 ENDDO
268
269 ! 63. traceurs
270
271 ! initialisation des tendances
272 dqfi=0.
273
274 DO iq=1, nqmx
275 iiq=niadv(iq)
276 DO l=1, llm
277 DO i=1, iim + 1
278 dqfi(i, 1, l, iiq) = d_qx(1, l, iq)
279 dqfi(i, jjm + 1, l, iiq) = d_qx(klon, l, iq)
280 ENDDO
281 DO j=2, jjm
282 ig0=1+(j-2)*iim
283 DO i=1, iim
284 dqfi(i, j, l, iiq) = d_qx(ig0+i, l, iq)
285 ENDDO
286 dqfi(iim + 1, j, l, iiq) = dqfi(1, j, l, iq)
287 ENDDO
288 ENDDO
289 ENDDO
290
291 ! 65. champ u:
292
293 DO l=1, llm
294 DO i=1, iim + 1
295 dufi(i, 1, l) = 0.
296 dufi(i, jjm + 1, l) = 0.
297 ENDDO
298
299 DO j=2, jjm
300 ig0=1+(j-2)*iim
301 DO i=1, iim-1
302 dufi(i, j, l)= 0.5*(d_u(ig0+i, l)+d_u(ig0+i+1, l))*cu_2d(i, j)
303 ENDDO
304 dufi(iim, j, l)= 0.5*(d_u(ig0+1, l)+d_u(ig0+iim, l))*cu_2d(iim, j)
305 dufi(iim + 1, j, l)=dufi(1, j, l)
306 ENDDO
307 ENDDO
308
309 ! 67. champ v:
310
311 DO l=1, llm
312 DO j=2, jjm-1
313 ig0=1+(j-2)*iim
314 DO i=1, iim
315 dvfi(i, j, l)= 0.5*(d_v(ig0+i, l)+d_v(ig0+i+iim, l))*cv_2d(i, j)
316 ENDDO
317 dvfi(iim + 1, j, l) = dvfi(1, j, l)
318 ENDDO
319 ENDDO
320
321 ! 68. champ v près des pôles:
322 ! v = U * cos(long) + V * SIN(long)
323
324 DO l=1, llm
325 DO i=1, iim
326 dvfi(i, 1, l)= d_u(1, l)*COS(rlonv(i))+d_v(1, l)*SIN(rlonv(i))
327 dvfi(i, jjm, l)=d_u(klon, l)*COS(rlonv(i)) +d_v(klon, l)*SIN(rlonv(i))
328 dvfi(i, 1, l)= 0.5*(dvfi(i, 1, l)+d_v(i+1, l))*cv_2d(i, 1)
329 dvfi(i, jjm, l)= 0.5 &
330 * (dvfi(i, jjm, l) + d_v(klon - iim - 1 + i, l)) * cv_2d(i, jjm)
331 ENDDO
332
333 dvfi(iim + 1, 1, l) = dvfi(1, 1, l)
334 dvfi(iim + 1, jjm, l)= dvfi(1, jjm, l)
335 ENDDO
336
337 END SUBROUTINE calfis
338
339 end module calfis_m

  ViewVC Help
Powered by ViewVC 1.1.21