/[lmdze]/trunk/phylmd/physiq.f
ViewVC logotype

Annotation of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 30 - (hide annotations)
Thu Apr 1 09:07:28 2010 UTC (14 years, 1 month ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 75385 byte(s)
Imported Source files of the external library "IOIPSL_Lionel" into
"libf/IOIPSL".

Split "cray.f90" into "scopy.f90" and "ssum.f90".

Rewrote "leapfrog" in order to have a clearer algorithmic structure.

1 guez 3 module physiq_m
2    
3     ! This module is clean: no C preprocessor directive, no include line.
4    
5     IMPLICIT none
6    
7     private
8     public physiq
9    
10     contains
11    
12 guez 12 SUBROUTINE physiq(nq, firstcal, lafin, rdayvrai, gmtime, pdtphys, paprs, &
13 guez 20 pplay, pphi, pphis, u, v, t, qx, omega, d_u, d_v, &
14 guez 3 d_t, d_qx, d_ps, dudyn, PVteta)
15    
16     ! From phylmd/physiq.F, v 1.22 2006/02/20 09:38:28
17    
18     ! Author : Z.X. Li (LMD/CNRS), date: 1993/08/18
19    
20     ! Objet: Moniteur general de la physique du modele
21     !AA Modifications quant aux traceurs :
22     !AA - uniformisation des parametrisations ds phytrac
23     !AA - stockage des moyennes des champs necessaires
24     !AA en mode traceur off-line
25    
26 guez 30 USE calendar, only: ymds2ju
27     USE histwrite_m, only: histwrite
28     USE histcom, only: histsync
29 guez 3 use dimens_m, only: jjm, iim, llm
30     use indicesol, only: nbsrf, is_ter, is_lic, is_sic, is_oce, &
31     clnsurf, epsfra
32     use dimphy, only: klon, nbtr
33 guez 23 use conf_gcm_m, only: raz_date, offline
34 guez 3 use dimsoil, only: nsoilmx
35 guez 23 use temps, only: itau_phy, day_ref, annee_ref
36 guez 15 use clesphys, only: ecrit_hf, ecrit_ins, ecrit_mth, &
37 guez 12 cdmmax, cdhmax, &
38     co2_ppm, ecrit_reg, ecrit_tra, ksta, ksta_ter, &
39     ok_kzmin
40     use clesphys2, only: iflag_con, ok_orolf, ok_orodr, nbapp_rad, &
41     cycle_diurne, new_oliq, soil_model
42     use iniprint, only: prt_level
43 guez 3 use abort_gcm_m, only: abort_gcm
44     use YOMCST, only: rcpd, rtt, rlvtt, rg, ra, rsigma, retv, romega
45     use comgeomphy
46     use ctherm
47     use phytrac_m, only: phytrac
48     use oasis_m
49     use radepsi
50     use radopt
51     use yoethf
52     use ini_hist, only: ini_histhf, ini_histday, ini_histins
53     use orbite_m, only: orbite, zenang
54     use phyetat0_m, only: phyetat0, rlat, rlon
55     use hgardfou_m, only: hgardfou
56     use conf_phys_m, only: conf_phys
57 guez 15 use phyredem_m, only: phyredem
58 guez 17 use qcheck_m, only: qcheck
59 guez 22 use ozonecm_m, only: ozonecm
60 guez 3
61     ! Declaration des constantes et des fonctions thermodynamiques :
62     use fcttre, only: thermcep, foeew, qsats, qsatl
63    
64     ! Variables argument:
65    
66 guez 17 INTEGER, intent(in):: nq ! nombre de traceurs (y compris vapeur d'eau)
67 guez 15
68     REAL, intent(in):: rdayvrai
69     ! (elapsed time since January 1st 0h of the starting year, in days)
70    
71 guez 3 REAL, intent(in):: gmtime ! heure de la journée en fraction de jour
72 guez 12 REAL, intent(in):: pdtphys ! pas d'integration pour la physique (seconde)
73 guez 7 LOGICAL, intent(in):: firstcal ! first call to "calfis"
74 guez 3 logical, intent(in):: lafin ! dernier passage
75    
76     REAL, intent(in):: paprs(klon, llm+1)
77     ! (pression pour chaque inter-couche, en Pa)
78 guez 15
79 guez 10 REAL, intent(in):: pplay(klon, llm)
80 guez 3 ! (input pression pour le mileu de chaque couche (en Pa))
81    
82     REAL pphi(klon, llm)
83     ! (input geopotentiel de chaque couche (g z) (reference sol))
84    
85     REAL pphis(klon) ! input geopotentiel du sol
86    
87     REAL u(klon, llm) ! input vitesse dans la direction X (de O a E) en m/s
88     REAL v(klon, llm) ! input vitesse Y (de S a N) en m/s
89     REAL t(klon, llm) ! input temperature (K)
90    
91 guez 12 REAL, intent(in):: qx(klon, llm, nq)
92     ! (humidite specifique (kg/kg) et fractions massiques des autres traceurs)
93 guez 3
94     REAL omega(klon, llm) ! input vitesse verticale en Pa/s
95     REAL d_u(klon, llm) ! output tendance physique de "u" (m/s/s)
96     REAL d_v(klon, llm) ! output tendance physique de "v" (m/s/s)
97     REAL d_t(klon, llm) ! output tendance physique de "t" (K/s)
98     REAL d_qx(klon, llm, nq) ! output tendance physique de "qx" (kg/kg/s)
99     REAL d_ps(klon) ! output tendance physique de la pression au sol
100    
101     INTEGER nbteta
102     PARAMETER(nbteta=3)
103    
104     REAL PVteta(klon, nbteta)
105     ! (output vorticite potentielle a des thetas constantes)
106    
107     LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
108     PARAMETER (ok_cvl=.TRUE.)
109     LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
110     PARAMETER (ok_gust=.FALSE.)
111    
112     LOGICAL check ! Verifier la conservation du modele en eau
113     PARAMETER (check=.FALSE.)
114     LOGICAL ok_stratus ! Ajouter artificiellement les stratus
115     PARAMETER (ok_stratus=.FALSE.)
116    
117     ! Parametres lies au coupleur OASIS:
118     INTEGER, SAVE :: npas, nexca
119     logical rnpb
120     parameter(rnpb=.true.)
121    
122 guez 12 character(len=6), save:: ocean
123     ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")
124    
125 guez 3 logical ok_ocean
126     SAVE ok_ocean
127    
128     !IM "slab" ocean
129     REAL tslab(klon) !Temperature du slab-ocean
130     SAVE tslab
131     REAL seaice(klon) !glace de mer (kg/m2)
132     SAVE seaice
133     REAL fluxo(klon) !flux turbulents ocean-glace de mer
134     REAL fluxg(klon) !flux turbulents ocean-atmosphere
135    
136     ! Modele thermique du sol, a activer pour le cycle diurne:
137 guez 17 logical, save:: ok_veget
138     LOGICAL, save:: ok_journe ! sortir le fichier journalier
139 guez 3
140     LOGICAL ok_mensuel ! sortir le fichier mensuel
141    
142     LOGICAL ok_instan ! sortir le fichier instantane
143     save ok_instan
144    
145     LOGICAL ok_region ! sortir le fichier regional
146     PARAMETER (ok_region=.FALSE.)
147    
148     ! pour phsystoke avec thermiques
149     REAL fm_therm(klon, llm+1)
150     REAL entr_therm(klon, llm)
151     real q2(klon, llm+1, nbsrf)
152     save q2
153    
154     INTEGER ivap ! indice de traceurs pour vapeur d'eau
155     PARAMETER (ivap=1)
156     INTEGER iliq ! indice de traceurs pour eau liquide
157     PARAMETER (iliq=2)
158    
159     REAL t_ancien(klon, llm), q_ancien(klon, llm)
160     SAVE t_ancien, q_ancien
161     LOGICAL ancien_ok
162     SAVE ancien_ok
163    
164     REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)
165     REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg/kg/s)
166    
167     real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
168    
169     !IM Amip2 PV a theta constante
170    
171     CHARACTER(LEN=3) ctetaSTD(nbteta)
172     DATA ctetaSTD/'350', '380', '405'/
173     REAL rtetaSTD(nbteta)
174     DATA rtetaSTD/350., 380., 405./
175    
176     !MI Amip2 PV a theta constante
177    
178     INTEGER klevp1
179     PARAMETER(klevp1=llm+1)
180    
181     REAL swdn0(klon, klevp1), swdn(klon, klevp1)
182     REAL swup0(klon, klevp1), swup(klon, klevp1)
183     SAVE swdn0, swdn, swup0, swup
184    
185     REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)
186     REAL lwup0(klon, klevp1), lwup(klon, klevp1)
187     SAVE lwdn0, lwdn, lwup0, lwup
188    
189     !IM Amip2
190     ! variables a une pression donnee
191    
192     integer nlevSTD
193     PARAMETER(nlevSTD=17)
194     real rlevSTD(nlevSTD)
195     DATA rlevSTD/100000., 92500., 85000., 70000., &
196     60000., 50000., 40000., 30000., 25000., 20000., &
197     15000., 10000., 7000., 5000., 3000., 2000., 1000./
198     CHARACTER(LEN=4) clevSTD(nlevSTD)
199     DATA clevSTD/'1000', '925 ', '850 ', '700 ', '600 ', &
200     '500 ', '400 ', '300 ', '250 ', '200 ', '150 ', '100 ', &
201     '70 ', '50 ', '30 ', '20 ', '10 '/
202    
203     ! prw: precipitable water
204     real prw(klon)
205    
206     ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
207     ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
208     REAL flwp(klon), fiwp(klon)
209     REAL flwc(klon, llm), fiwc(klon, llm)
210    
211 guez 23 INTEGER kmax, lmax
212 guez 3 PARAMETER(kmax=8, lmax=8)
213     INTEGER kmaxm1, lmaxm1
214     PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
215    
216     REAL zx_tau(kmaxm1), zx_pc(lmaxm1)
217     DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
218     DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
219    
220     ! cldtopres pression au sommet des nuages
221     REAL cldtopres(lmaxm1)
222     DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
223    
224     ! taulev: numero du niveau de tau dans les sorties ISCCP
225     CHARACTER(LEN=4) taulev(kmaxm1)
226    
227     DATA taulev/'tau0', 'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'/
228     CHARACTER(LEN=3) pclev(lmaxm1)
229     DATA pclev/'pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7'/
230    
231     CHARACTER(LEN=28) cnameisccp(lmaxm1, kmaxm1)
232     DATA cnameisccp/'pc< 50hPa, tau< 0.3', 'pc= 50-180hPa, tau< 0.3', &
233     'pc= 180-310hPa, tau< 0.3', 'pc= 310-440hPa, tau< 0.3', &
234     'pc= 440-560hPa, tau< 0.3', 'pc= 560-680hPa, tau< 0.3', &
235     'pc= 680-800hPa, tau< 0.3', 'pc< 50hPa, tau= 0.3-1.3', &
236     'pc= 50-180hPa, tau= 0.3-1.3', 'pc= 180-310hPa, tau= 0.3-1.3', &
237     'pc= 310-440hPa, tau= 0.3-1.3', 'pc= 440-560hPa, tau= 0.3-1.3', &
238     'pc= 560-680hPa, tau= 0.3-1.3', 'pc= 680-800hPa, tau= 0.3-1.3', &
239     'pc< 50hPa, tau= 1.3-3.6', 'pc= 50-180hPa, tau= 1.3-3.6', &
240     'pc= 180-310hPa, tau= 1.3-3.6', 'pc= 310-440hPa, tau= 1.3-3.6', &
241     'pc= 440-560hPa, tau= 1.3-3.6', 'pc= 560-680hPa, tau= 1.3-3.6', &
242     'pc= 680-800hPa, tau= 1.3-3.6', 'pc< 50hPa, tau= 3.6-9.4', &
243     'pc= 50-180hPa, tau= 3.6-9.4', 'pc= 180-310hPa, tau= 3.6-9.4', &
244     'pc= 310-440hPa, tau= 3.6-9.4', 'pc= 440-560hPa, tau= 3.6-9.4', &
245     'pc= 560-680hPa, tau= 3.6-9.4', 'pc= 680-800hPa, tau= 3.6-9.4', &
246     'pc< 50hPa, tau= 9.4-23', 'pc= 50-180hPa, tau= 9.4-23', &
247     'pc= 180-310hPa, tau= 9.4-23', 'pc= 310-440hPa, tau= 9.4-23', &
248     'pc= 440-560hPa, tau= 9.4-23', 'pc= 560-680hPa, tau= 9.4-23', &
249     'pc= 680-800hPa, tau= 9.4-23', 'pc< 50hPa, tau= 23-60', &
250     'pc= 50-180hPa, tau= 23-60', 'pc= 180-310hPa, tau= 23-60', &
251     'pc= 310-440hPa, tau= 23-60', 'pc= 440-560hPa, tau= 23-60', &
252     'pc= 560-680hPa, tau= 23-60', 'pc= 680-800hPa, tau= 23-60', &
253     'pc< 50hPa, tau> 60.', 'pc= 50-180hPa, tau> 60.', &
254     'pc= 180-310hPa, tau> 60.', 'pc= 310-440hPa, tau> 60.', &
255     'pc= 440-560hPa, tau> 60.', 'pc= 560-680hPa, tau> 60.', &
256     'pc= 680-800hPa, tau> 60.'/
257    
258     !IM ISCCP simulator v3.4
259    
260     integer nid_hf, nid_hf3d
261     save nid_hf, nid_hf3d
262    
263     ! Variables propres a la physique
264    
265     INTEGER, save:: radpas
266     ! (Radiative transfer computations are made every "radpas" call to
267     ! "physiq".)
268    
269     REAL radsol(klon)
270     SAVE radsol ! bilan radiatif au sol calcule par code radiatif
271    
272 guez 7 INTEGER, SAVE:: itap ! number of calls to "physiq"
273 guez 3
274     REAL ftsol(klon, nbsrf)
275     SAVE ftsol ! temperature du sol
276    
277     REAL ftsoil(klon, nsoilmx, nbsrf)
278     SAVE ftsoil ! temperature dans le sol
279    
280     REAL fevap(klon, nbsrf)
281     SAVE fevap ! evaporation
282     REAL fluxlat(klon, nbsrf)
283     SAVE fluxlat
284    
285     REAL fqsurf(klon, nbsrf)
286     SAVE fqsurf ! humidite de l'air au contact de la surface
287    
288     REAL qsol(klon)
289     SAVE qsol ! hauteur d'eau dans le sol
290    
291     REAL fsnow(klon, nbsrf)
292     SAVE fsnow ! epaisseur neigeuse
293    
294     REAL falbe(klon, nbsrf)
295     SAVE falbe ! albedo par type de surface
296     REAL falblw(klon, nbsrf)
297     SAVE falblw ! albedo par type de surface
298    
299 guez 13 ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :
300     REAL, save:: zmea(klon) ! orographie moyenne
301     REAL, save:: zstd(klon) ! deviation standard de l'OESM
302     REAL, save:: zsig(klon) ! pente de l'OESM
303     REAL, save:: zgam(klon) ! anisotropie de l'OESM
304     REAL, save:: zthe(klon) ! orientation de l'OESM
305     REAL, save:: zpic(klon) ! Maximum de l'OESM
306     REAL, save:: zval(klon) ! Minimum de l'OESM
307     REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
308 guez 3
309     REAL zulow(klon), zvlow(klon)
310    
311     INTEGER igwd, idx(klon), itest(klon)
312    
313     REAL agesno(klon, nbsrf)
314     SAVE agesno ! age de la neige
315    
316     REAL run_off_lic_0(klon)
317     SAVE run_off_lic_0
318     !KE43
319     ! Variables liees a la convection de K. Emanuel (sb):
320    
321     REAL bas, top ! cloud base and top levels
322     SAVE bas
323     SAVE top
324    
325     REAL Ma(klon, llm) ! undilute upward mass flux
326     SAVE Ma
327     REAL qcondc(klon, llm) ! in-cld water content from convect
328     SAVE qcondc
329     REAL ema_work1(klon, llm), ema_work2(klon, llm)
330     SAVE ema_work1, ema_work2
331    
332     REAL wd(klon) ! sb
333     SAVE wd ! sb
334    
335     ! Variables locales pour la couche limite (al1):
336    
337     ! Variables locales:
338    
339     REAL cdragh(klon) ! drag coefficient pour T and Q
340     REAL cdragm(klon) ! drag coefficient pour vent
341    
342     !AA Pour phytrac
343     REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
344     REAL yu1(klon) ! vents dans la premiere couche U
345     REAL yv1(klon) ! vents dans la premiere couche V
346     REAL ffonte(klon, nbsrf) !Flux thermique utilise pour fondre la neige
347     REAL fqcalving(klon, nbsrf) !Flux d'eau "perdue" par la surface
348     ! !et necessaire pour limiter la
349     ! !hauteur de neige, en kg/m2/s
350     REAL zxffonte(klon), zxfqcalving(klon)
351    
352     REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
353     save pfrac_impa
354     REAL pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation
355     save pfrac_nucl
356     REAL pfrac_1nucl(klon, llm)! Produits des coefs lessi nucl (alpha = 1)
357     save pfrac_1nucl
358     REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
359     REAL frac_nucl(klon, llm) ! idem (nucleation)
360    
361     !AA
362     REAL rain_fall(klon) ! pluie
363     REAL snow_fall(klon) ! neige
364     save snow_fall, rain_fall
365     !IM cf FH pour Tiedtke 080604
366     REAL rain_tiedtke(klon), snow_tiedtke(klon)
367    
368     REAL evap(klon), devap(klon) ! evaporation et sa derivee
369     REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
370     REAL dlw(klon) ! derivee infra rouge
371     SAVE dlw
372     REAL bils(klon) ! bilan de chaleur au sol
373     REAL fder(klon) ! Derive de flux (sensible et latente)
374     save fder
375     REAL ve(klon) ! integr. verticale du transport meri. de l'energie
376     REAL vq(klon) ! integr. verticale du transport meri. de l'eau
377     REAL ue(klon) ! integr. verticale du transport zonal de l'energie
378     REAL uq(klon) ! integr. verticale du transport zonal de l'eau
379    
380     REAL frugs(klon, nbsrf) ! longueur de rugosite
381     save frugs
382     REAL zxrugs(klon) ! longueur de rugosite
383    
384     ! Conditions aux limites
385    
386     INTEGER julien
387    
388 guez 7 INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
389 guez 3 REAL pctsrf(klon, nbsrf)
390     !IM
391     REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE
392    
393     SAVE pctsrf ! sous-fraction du sol
394     REAL albsol(klon)
395     SAVE albsol ! albedo du sol total
396     REAL albsollw(klon)
397     SAVE albsollw ! albedo du sol total
398    
399 guez 17 REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
400 guez 3
401     ! Declaration des procedures appelees
402    
403     EXTERNAL alboc ! calculer l'albedo sur ocean
404     EXTERNAL ajsec ! ajustement sec
405     EXTERNAL clmain ! couche limite
406     !KE43
407     EXTERNAL conema3 ! convect4.3
408     EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
409     EXTERNAL nuage ! calculer les proprietes radiatives
410     EXTERNAL radlwsw ! rayonnements solaire et infrarouge
411     EXTERNAL transp ! transport total de l'eau et de l'energie
412    
413     ! Variables locales
414    
415     real clwcon(klon, llm), rnebcon(klon, llm)
416     real clwcon0(klon, llm), rnebcon0(klon, llm)
417    
418     save rnebcon, clwcon
419    
420     REAL rhcl(klon, llm) ! humiditi relative ciel clair
421     REAL dialiq(klon, llm) ! eau liquide nuageuse
422     REAL diafra(klon, llm) ! fraction nuageuse
423     REAL cldliq(klon, llm) ! eau liquide nuageuse
424     REAL cldfra(klon, llm) ! fraction nuageuse
425     REAL cldtau(klon, llm) ! epaisseur optique
426     REAL cldemi(klon, llm) ! emissivite infrarouge
427    
428     REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite
429     REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur
430     REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u
431     REAL fluxv(klon, llm, nbsrf) ! flux turbulent de vitesse v
432    
433     REAL zxfluxt(klon, llm)
434     REAL zxfluxq(klon, llm)
435     REAL zxfluxu(klon, llm)
436     REAL zxfluxv(klon, llm)
437    
438     REAL heat(klon, llm) ! chauffage solaire
439     REAL heat0(klon, llm) ! chauffage solaire ciel clair
440     REAL cool(klon, llm) ! refroidissement infrarouge
441     REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair
442     REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
443     real sollwdown(klon) ! downward LW flux at surface
444     REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
445     REAL albpla(klon)
446     REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
447     REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
448     ! Le rayonnement n'est pas calcule tous les pas, il faut donc
449     ! sauvegarder les sorties du rayonnement
450     SAVE heat, cool, albpla, topsw, toplw, solsw, sollw, sollwdown
451     SAVE topsw0, toplw0, solsw0, sollw0, heat0, cool0
452    
453     INTEGER itaprad
454     SAVE itaprad
455    
456     REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
457     REAL conv_t(klon, llm) ! convergence de la temperature(K/s)
458    
459     REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut
460     REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree
461    
462     REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
463    
464     REAL dist, rmu0(klon), fract(klon)
465     REAL zdtime ! pas de temps du rayonnement (s)
466     real zlongi
467    
468     REAL z_avant(klon), z_apres(klon), z_factor(klon)
469     LOGICAL zx_ajustq
470    
471     REAL za, zb
472     REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp
473     real zqsat(klon, llm)
474     INTEGER i, k, iq, nsrf
475     REAL t_coup
476     PARAMETER (t_coup=234.0)
477    
478     REAL zphi(klon, llm)
479    
480     !IM cf. AM Variables locales pour la CLA (hbtm2)
481    
482     REAL pblh(klon, nbsrf) ! Hauteur de couche limite
483     REAL plcl(klon, nbsrf) ! Niveau de condensation de la CLA
484     REAL capCL(klon, nbsrf) ! CAPE de couche limite
485     REAL oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
486     REAL cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
487     REAL pblt(klon, nbsrf) ! T a la Hauteur de couche limite
488     REAL therm(klon, nbsrf)
489     REAL trmb1(klon, nbsrf) ! deep_cape
490     REAL trmb2(klon, nbsrf) ! inhibition
491     REAL trmb3(klon, nbsrf) ! Point Omega
492     ! Grdeurs de sorties
493     REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
494     REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
495     REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
496     REAL s_trmb3(klon)
497    
498     ! Variables locales pour la convection de K. Emanuel (sb):
499    
500     REAL upwd(klon, llm) ! saturated updraft mass flux
501     REAL dnwd(klon, llm) ! saturated downdraft mass flux
502     REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
503     REAL tvp(klon, llm) ! virtual temp of lifted parcel
504     REAL cape(klon) ! CAPE
505     SAVE cape
506    
507     REAL pbase(klon) ! cloud base pressure
508     SAVE pbase
509     REAL bbase(klon) ! cloud base buoyancy
510     SAVE bbase
511     REAL rflag(klon) ! flag fonctionnement de convect
512     INTEGER iflagctrl(klon) ! flag fonctionnement de convect
513     ! -- convect43:
514     INTEGER ntra ! nb traceurs pour convect4.3
515     REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)
516     REAL dplcldt(klon), dplcldr(klon)
517    
518     ! Variables du changement
519    
520     ! con: convection
521     ! lsc: condensation a grande echelle (Large-Scale-Condensation)
522     ! ajs: ajustement sec
523     ! eva: evaporation de l'eau liquide nuageuse
524     ! vdf: couche limite (Vertical DiFfusion)
525     REAL d_t_con(klon, llm), d_q_con(klon, llm)
526     REAL d_u_con(klon, llm), d_v_con(klon, llm)
527     REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
528     REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
529     REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
530     REAL rneb(klon, llm)
531    
532     REAL pmfu(klon, llm), pmfd(klon, llm)
533     REAL pen_u(klon, llm), pen_d(klon, llm)
534     REAL pde_u(klon, llm), pde_d(klon, llm)
535     INTEGER kcbot(klon), kctop(klon), kdtop(klon)
536     REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1)
537     REAL prfl(klon, llm+1), psfl(klon, llm+1)
538    
539     INTEGER ibas_con(klon), itop_con(klon)
540    
541     SAVE ibas_con, itop_con
542    
543     REAL rain_con(klon), rain_lsc(klon)
544     REAL snow_con(klon), snow_lsc(klon)
545     REAL d_ts(klon, nbsrf)
546    
547     REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
548     REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
549    
550     REAL d_u_oro(klon, llm), d_v_oro(klon, llm)
551     REAL d_t_oro(klon, llm)
552     REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
553     REAL d_t_lif(klon, llm)
554    
555     REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm)
556     real ratqsbas, ratqshaut
557     save ratqsbas, ratqshaut, ratqs
558    
559     ! Parametres lies au nouveau schema de nuages (SB, PDF)
560 guez 17 real, save:: fact_cldcon
561     real, save:: facttemps
562 guez 3 logical ok_newmicro
563     save ok_newmicro
564     real facteur
565    
566     integer iflag_cldcon
567     save iflag_cldcon
568    
569     logical ptconv(klon, llm)
570    
571     ! Variables locales pour effectuer les appels en serie
572    
573     REAL t_seri(klon, llm), q_seri(klon, llm)
574     REAL ql_seri(klon, llm), qs_seri(klon, llm)
575     REAL u_seri(klon, llm), v_seri(klon, llm)
576    
577     REAL tr_seri(klon, llm, nbtr)
578     REAL d_tr(klon, llm, nbtr)
579    
580     REAL zx_rh(klon, llm)
581    
582     REAL zustrdr(klon), zvstrdr(klon)
583     REAL zustrli(klon), zvstrli(klon)
584     REAL zustrph(klon), zvstrph(klon)
585     REAL aam, torsfc
586    
587     REAL dudyn(iim+1, jjm + 1, llm)
588    
589     REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
590     REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
591    
592 guez 15 INTEGER, SAVE:: nid_day, nid_ins
593 guez 3
594     REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
595     REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
596     REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
597     REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
598    
599     REAL zsto
600    
601     character(len=20) modname
602     character(len=80) abort_message
603     logical ok_sync
604     real date0
605    
606     ! Variables liees au bilan d'energie et d'enthalpi
607     REAL ztsol(klon)
608     REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
609     REAL d_h_vcol_phy
610     REAL fs_bound, fq_bound
611     SAVE d_h_vcol_phy
612     REAL zero_v(klon)
613     CHARACTER(LEN=15) ztit
614     INTEGER ip_ebil ! PRINT level for energy conserv. diag.
615     SAVE ip_ebil
616     DATA ip_ebil/0/
617 guez 18 INTEGER, SAVE:: if_ebil ! level for energy conservation diagnostics
618 guez 3 !+jld ec_conser
619     REAL d_t_ec(klon, llm) ! tendance du a la conersion Ec -> E thermique
620     REAL ZRCPD
621     !-jld ec_conser
622     !IM: t2m, q2m, u10m, v10m
623     REAL t2m(klon, nbsrf), q2m(klon, nbsrf) !temperature, humidite a 2m
624     REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m
625     REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
626     REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
627     !jq Aerosol effects (Johannes Quaas, 27/11/2003)
628     REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3]
629    
630     REAL sulfate_pi(klon, llm)
631     ! (SO4 aerosol concentration [ug/m3] (pre-industrial value))
632     SAVE sulfate_pi
633    
634     REAL cldtaupi(klon, llm)
635     ! (Cloud optical thickness for pre-industrial (pi) aerosols)
636    
637     REAL re(klon, llm) ! Cloud droplet effective radius
638     REAL fl(klon, llm) ! denominator of re
639    
640     ! Aerosol optical properties
641     REAL tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
642     REAL cg_ae(klon, llm, 2)
643    
644     REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
645     ! ok_ade=T -ADE=topswad-topsw
646    
647     REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
648     ! ok_aie=T ->
649     ! ok_ade=T -AIE=topswai-topswad
650     ! ok_ade=F -AIE=topswai-topsw
651    
652     REAL aerindex(klon) ! POLDER aerosol index
653    
654     ! Parameters
655     LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
656     REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
657    
658     SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
659     SAVE u10m
660     SAVE v10m
661     SAVE t2m
662     SAVE q2m
663     SAVE ffonte
664     SAVE fqcalving
665     SAVE piz_ae
666     SAVE tau_ae
667     SAVE cg_ae
668     SAVE rain_con
669     SAVE snow_con
670     SAVE topswai
671     SAVE topswad
672     SAVE solswai
673     SAVE solswad
674     SAVE d_u_con
675     SAVE d_v_con
676     SAVE rnebcon0
677     SAVE clwcon0
678     SAVE pblh
679     SAVE plcl
680     SAVE capCL
681     SAVE oliqCL
682     SAVE cteiCL
683     SAVE pblt
684     SAVE therm
685     SAVE trmb1
686     SAVE trmb2
687     SAVE trmb3
688    
689 guez 17 real zmasse(klon, llm)
690     ! (column-density of mass of air in a cell, in kg m-2)
691    
692     real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
693    
694 guez 3 !----------------------------------------------------------------
695    
696     modname = 'physiq'
697     IF (if_ebil >= 1) THEN
698     DO i=1, klon
699     zero_v(i)=0.
700     END DO
701     END IF
702     ok_sync=.TRUE.
703 guez 7 IF (nq < 2) THEN
704 guez 3 abort_message = 'eaux vapeur et liquide sont indispensables'
705 guez 13 CALL abort_gcm(modname, abort_message, 1)
706 guez 3 ENDIF
707    
708 guez 7 test_firstcal: IF (firstcal) THEN
709 guez 3 ! initialiser
710 guez 13 u10m=0.
711     v10m=0.
712     t2m=0.
713     q2m=0.
714     ffonte=0.
715     fqcalving=0.
716 guez 3 piz_ae(:, :, :)=0.
717     tau_ae(:, :, :)=0.
718     cg_ae(:, :, :)=0.
719     rain_con(:)=0.
720     snow_con(:)=0.
721     bl95_b0=0.
722     bl95_b1=0.
723     topswai(:)=0.
724     topswad(:)=0.
725     solswai(:)=0.
726     solswad(:)=0.
727    
728 guez 13 d_u_con = 0.0
729     d_v_con = 0.0
730     rnebcon0 = 0.0
731     clwcon0 = 0.0
732     rnebcon = 0.0
733     clwcon = 0.0
734 guez 3
735 guez 13 pblh =0. ! Hauteur de couche limite
736     plcl =0. ! Niveau de condensation de la CLA
737     capCL =0. ! CAPE de couche limite
738     oliqCL =0. ! eau_liqu integree de couche limite
739     cteiCL =0. ! cloud top instab. crit. couche limite
740     pblt =0. ! T a la Hauteur de couche limite
741     therm =0.
742     trmb1 =0. ! deep_cape
743     trmb2 =0. ! inhibition
744     trmb3 =0. ! Point Omega
745 guez 3
746     IF (if_ebil >= 1) d_h_vcol_phy=0.
747    
748     ! appel a la lecture du run.def physique
749    
750     call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, &
751     ok_instan, fact_cldcon, facttemps, ok_newmicro, &
752     iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &
753     ok_ade, ok_aie, &
754     bl95_b0, bl95_b1, &
755     iflag_thermals, nsplit_thermals)
756    
757     ! Initialiser les compteurs:
758    
759     frugs = 0.
760     itap = 0
761     itaprad = 0
762 guez 12 CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &
763     seaice, fqsurf, qsol, fsnow, &
764 guez 3 falbe, falblw, fevap, rain_fall, snow_fall, solsw, sollwdown, &
765 guez 13 dlw, radsol, frugs, agesno, &
766     zmea, zstd, zsig, zgam, zthe, zpic, zval, &
767 guez 3 t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon, &
768     run_off_lic_0)
769    
770     ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
771     q2(:, :, :)=1.e-8
772    
773 guez 12 radpas = NINT( 86400. / pdtphys / nbapp_rad)
774 guez 3
775     ! on remet le calendrier a zero
776 guez 15 IF (raz_date) itau_phy = 0
777 guez 3
778 guez 12 PRINT *, 'cycle_diurne = ', cycle_diurne
779 guez 3
780     IF(ocean.NE.'force ') THEN
781     ok_ocean=.TRUE.
782     ENDIF
783    
784 guez 12 CALL printflag(radpas, ok_ocean, ok_oasis, ok_journe, ok_instan, &
785     ok_region)
786 guez 3
787 guez 12 IF (pdtphys*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
788     print *,'Nbre d appels au rayonnement insuffisant'
789     print *,"Au minimum 4 appels par jour si cycle diurne"
790 guez 3 abort_message='Nbre d appels au rayonnement insuffisant'
791     call abort_gcm(modname, abort_message, 1)
792     ENDIF
793 guez 12 print *,"Clef pour la convection, iflag_con=", iflag_con
794     print *,"Clef pour le driver de la convection, ok_cvl=", &
795 guez 3 ok_cvl
796    
797     ! Initialisation pour la convection de K.E. (sb):
798     IF (iflag_con >= 3) THEN
799    
800 guez 12 print *,"*** Convection de Kerry Emanuel 4.3 "
801 guez 3
802     !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>BEG
803     DO i = 1, klon
804     ibas_con(i) = 1
805     itop_con(i) = 1
806     ENDDO
807     !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>END
808    
809     ENDIF
810    
811     IF (ok_orodr) THEN
812 guez 13 rugoro = MAX(1e-5, zstd * zsig / 2)
813 guez 3 CALL SUGWD(klon, llm, paprs, pplay)
814 guez 13 else
815     rugoro = 0.
816 guez 3 ENDIF
817    
818 guez 12 lmt_pas = NINT(86400. / pdtphys) ! tous les jours
819 guez 7 print *, 'Number of time steps of "physics" per day: ', lmt_pas
820 guez 3
821 guez 12 ecrit_ins = NINT(ecrit_ins/pdtphys)
822     ecrit_hf = NINT(ecrit_hf/pdtphys)
823     ecrit_mth = NINT(ecrit_mth/pdtphys)
824     ecrit_tra = NINT(86400.*ecrit_tra/pdtphys)
825     ecrit_reg = NINT(ecrit_reg/pdtphys)
826 guez 3
827     ! Initialiser le couplage si necessaire
828    
829     npas = 0
830     nexca = 0
831    
832 guez 12 print *,'AVANT HIST IFLAG_CON=', iflag_con
833 guez 3
834     ! Initialisation des sorties
835    
836 guez 20 call ini_histhf(pdtphys, nid_hf, nid_hf3d)
837     call ini_histday(pdtphys, ok_journe, nid_day, nq)
838     call ini_histins(pdtphys, ok_instan, nid_ins)
839 guez 3 CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)
840     !XXXPB Positionner date0 pour initialisation de ORCHIDEE
841     WRITE(*, *) 'physiq date0 : ', date0
842 guez 7 ENDIF test_firstcal
843 guez 3
844     ! Mettre a zero des variables de sortie (pour securite)
845    
846     DO i = 1, klon
847     d_ps(i) = 0.0
848     ENDDO
849     DO k = 1, llm
850     DO i = 1, klon
851     d_t(i, k) = 0.0
852     d_u(i, k) = 0.0
853     d_v(i, k) = 0.0
854     ENDDO
855     ENDDO
856     DO iq = 1, nq
857     DO k = 1, llm
858     DO i = 1, klon
859     d_qx(i, k, iq) = 0.0
860     ENDDO
861     ENDDO
862     ENDDO
863 guez 13 da=0.
864     mp=0.
865 guez 3 phi(:, :, :)=0.
866    
867     ! Ne pas affecter les valeurs entrees de u, v, h, et q
868    
869     DO k = 1, llm
870     DO i = 1, klon
871     t_seri(i, k) = t(i, k)
872     u_seri(i, k) = u(i, k)
873     v_seri(i, k) = v(i, k)
874     q_seri(i, k) = qx(i, k, ivap)
875     ql_seri(i, k) = qx(i, k, iliq)
876     qs_seri(i, k) = 0.
877     ENDDO
878     ENDDO
879     IF (nq >= 3) THEN
880     tr_seri(:, :, :nq-2) = qx(:, :, 3:nq)
881     ELSE
882     tr_seri(:, :, 1) = 0.
883     ENDIF
884    
885     DO i = 1, klon
886     ztsol(i) = 0.
887     ENDDO
888     DO nsrf = 1, nbsrf
889     DO i = 1, klon
890     ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
891     ENDDO
892     ENDDO
893    
894     IF (if_ebil >= 1) THEN
895     ztit='after dynamic'
896 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, pdtphys &
897 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
898 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
899     ! Comme les tendances de la physique sont ajoute dans la dynamique,
900     ! on devrait avoir que la variation d'entalpie par la dynamique
901     ! est egale a la variation de la physique au pas de temps precedent.
902     ! Donc la somme de ces 2 variations devrait etre nulle.
903     call diagphy(airephy, ztit, ip_ebil &
904     , zero_v, zero_v, zero_v, zero_v, zero_v &
905     , zero_v, zero_v, zero_v, ztsol &
906     , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
907     , fs_bound, fq_bound )
908     END IF
909    
910     ! Diagnostiquer la tendance dynamique
911    
912     IF (ancien_ok) THEN
913     DO k = 1, llm
914     DO i = 1, klon
915 guez 12 d_t_dyn(i, k) = (t_seri(i, k)-t_ancien(i, k))/pdtphys
916     d_q_dyn(i, k) = (q_seri(i, k)-q_ancien(i, k))/pdtphys
917 guez 3 ENDDO
918     ENDDO
919     ELSE
920     DO k = 1, llm
921     DO i = 1, klon
922     d_t_dyn(i, k) = 0.0
923     d_q_dyn(i, k) = 0.0
924     ENDDO
925     ENDDO
926     ancien_ok = .TRUE.
927     ENDIF
928    
929     ! Ajouter le geopotentiel du sol:
930    
931     DO k = 1, llm
932     DO i = 1, klon
933     zphi(i, k) = pphi(i, k) + pphis(i)
934     ENDDO
935     ENDDO
936    
937     ! Verifier les temperatures
938    
939     CALL hgardfou(t_seri, ftsol)
940    
941     ! Incrementer le compteur de la physique
942    
943 guez 7 itap = itap + 1
944     julien = MOD(NINT(rdayvrai), 360)
945 guez 3 if (julien == 0) julien = 360
946    
947 guez 17 forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
948    
949 guez 3 ! Mettre en action les conditions aux limites (albedo, sst, etc.).
950     ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
951    
952 guez 18 if (nq >= 5) then
953     wo = qx(:, :, 5) * zmasse / dobson_u / 1e3
954     else IF (MOD(itap - 1, lmt_pas) == 0) THEN
955 guez 22 wo = ozonecm(REAL(julien), paprs)
956 guez 3 ENDIF
957    
958     ! Re-evaporer l'eau liquide nuageuse
959    
960     DO k = 1, llm ! re-evaporation de l'eau liquide nuageuse
961     DO i = 1, klon
962     zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i, k))
963     zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i, k))
964     zdelta = MAX(0., SIGN(1., RTT-t_seri(i, k)))
965     zb = MAX(0.0, ql_seri(i, k))
966     za = - MAX(0.0, ql_seri(i, k)) &
967     * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
968     t_seri(i, k) = t_seri(i, k) + za
969     q_seri(i, k) = q_seri(i, k) + zb
970     ql_seri(i, k) = 0.0
971     ENDDO
972     ENDDO
973    
974     IF (if_ebil >= 2) THEN
975     ztit='after reevap'
976 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 1, pdtphys &
977 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
978 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
979     call diagphy(airephy, ztit, ip_ebil &
980     , zero_v, zero_v, zero_v, zero_v, zero_v &
981     , zero_v, zero_v, zero_v, ztsol &
982     , d_h_vcol, d_qt, d_ec &
983     , fs_bound, fq_bound )
984    
985     END IF
986    
987     ! Appeler la diffusion verticale (programme de couche limite)
988    
989     DO i = 1, klon
990     zxrugs(i) = 0.0
991     ENDDO
992     DO nsrf = 1, nbsrf
993     DO i = 1, klon
994     frugs(i, nsrf) = MAX(frugs(i, nsrf), 0.000015)
995     ENDDO
996     ENDDO
997     DO nsrf = 1, nbsrf
998     DO i = 1, klon
999     zxrugs(i) = zxrugs(i) + frugs(i, nsrf)*pctsrf(i, nsrf)
1000     ENDDO
1001     ENDDO
1002    
1003     ! calculs necessaires au calcul de l'albedo dans l'interface
1004    
1005     CALL orbite(REAL(julien), zlongi, dist)
1006     IF (cycle_diurne) THEN
1007 guez 12 zdtime = pdtphys * REAL(radpas)
1008 guez 3 CALL zenang(zlongi, gmtime, zdtime, rmu0, fract)
1009     ELSE
1010     rmu0 = -999.999
1011     ENDIF
1012    
1013     ! Calcul de l'abedo moyen par maille
1014     albsol(:)=0.
1015     albsollw(:)=0.
1016     DO nsrf = 1, nbsrf
1017     DO i = 1, klon
1018     albsol(i) = albsol(i) + falbe(i, nsrf) * pctsrf(i, nsrf)
1019     albsollw(i) = albsollw(i) + falblw(i, nsrf) * pctsrf(i, nsrf)
1020     ENDDO
1021     ENDDO
1022    
1023     ! Repartition sous maille des flux LW et SW
1024     ! Repartition du longwave par sous-surface linearisee
1025    
1026     DO nsrf = 1, nbsrf
1027     DO i = 1, klon
1028     fsollw(i, nsrf) = sollw(i) &
1029     + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))
1030     fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))
1031     ENDDO
1032     ENDDO
1033    
1034     fder = dlw
1035    
1036 guez 12 CALL clmain(pdtphys, itap, date0, pctsrf, pctsrf_new, &
1037 guez 3 t_seri, q_seri, u_seri, v_seri, &
1038     julien, rmu0, co2_ppm, &
1039     ok_veget, ocean, npas, nexca, ftsol, &
1040     soil_model, cdmmax, cdhmax, &
1041     ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &
1042     paprs, pplay, fsnow, fqsurf, fevap, falbe, falblw, &
1043     fluxlat, rain_fall, snow_fall, &
1044     fsolsw, fsollw, sollwdown, fder, &
1045     rlon, rlat, cuphy, cvphy, frugs, &
1046 guez 7 firstcal, lafin, agesno, rugoro, &
1047 guez 3 d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, &
1048     fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, &
1049     q2, dsens, devap, &
1050     ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &
1051     pblh, capCL, oliqCL, cteiCL, pblT, &
1052     therm, trmb1, trmb2, trmb3, plcl, &
1053     fqcalving, ffonte, run_off_lic_0, &
1054     fluxo, fluxg, tslab, seaice)
1055    
1056     !XXX Incrementation des flux
1057    
1058     zxfluxt=0.
1059     zxfluxq=0.
1060     zxfluxu=0.
1061     zxfluxv=0.
1062     DO nsrf = 1, nbsrf
1063     DO k = 1, llm
1064     DO i = 1, klon
1065     zxfluxt(i, k) = zxfluxt(i, k) + &
1066     fluxt(i, k, nsrf) * pctsrf( i, nsrf)
1067     zxfluxq(i, k) = zxfluxq(i, k) + &
1068     fluxq(i, k, nsrf) * pctsrf( i, nsrf)
1069     zxfluxu(i, k) = zxfluxu(i, k) + &
1070     fluxu(i, k, nsrf) * pctsrf( i, nsrf)
1071     zxfluxv(i, k) = zxfluxv(i, k) + &
1072     fluxv(i, k, nsrf) * pctsrf( i, nsrf)
1073     END DO
1074     END DO
1075     END DO
1076     DO i = 1, klon
1077     sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
1078     evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol
1079     fder(i) = dlw(i) + dsens(i) + devap(i)
1080     ENDDO
1081    
1082     DO k = 1, llm
1083     DO i = 1, klon
1084     t_seri(i, k) = t_seri(i, k) + d_t_vdf(i, k)
1085     q_seri(i, k) = q_seri(i, k) + d_q_vdf(i, k)
1086     u_seri(i, k) = u_seri(i, k) + d_u_vdf(i, k)
1087     v_seri(i, k) = v_seri(i, k) + d_v_vdf(i, k)
1088     ENDDO
1089     ENDDO
1090    
1091     IF (if_ebil >= 2) THEN
1092     ztit='after clmain'
1093 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1094 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1095 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1096     call diagphy(airephy, ztit, ip_ebil &
1097     , zero_v, zero_v, zero_v, zero_v, sens &
1098     , evap, zero_v, zero_v, ztsol &
1099     , d_h_vcol, d_qt, d_ec &
1100     , fs_bound, fq_bound )
1101     END IF
1102    
1103     ! Incrementer la temperature du sol
1104    
1105     DO i = 1, klon
1106     zxtsol(i) = 0.0
1107     zxfluxlat(i) = 0.0
1108    
1109     zt2m(i) = 0.0
1110     zq2m(i) = 0.0
1111     zu10m(i) = 0.0
1112     zv10m(i) = 0.0
1113     zxffonte(i) = 0.0
1114     zxfqcalving(i) = 0.0
1115    
1116     s_pblh(i) = 0.0
1117     s_lcl(i) = 0.0
1118     s_capCL(i) = 0.0
1119     s_oliqCL(i) = 0.0
1120     s_cteiCL(i) = 0.0
1121     s_pblT(i) = 0.0
1122     s_therm(i) = 0.0
1123     s_trmb1(i) = 0.0
1124     s_trmb2(i) = 0.0
1125     s_trmb3(i) = 0.0
1126    
1127     IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) + &
1128     pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA) &
1129     THEN
1130     WRITE(*, *) 'physiq : pb sous surface au point ', i, &
1131     pctsrf(i, 1 : nbsrf)
1132     ENDIF
1133     ENDDO
1134     DO nsrf = 1, nbsrf
1135     DO i = 1, klon
1136     ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)
1137     zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
1138     zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)
1139    
1140     zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)
1141     zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)
1142     zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)
1143     zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)
1144     zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)
1145     zxfqcalving(i) = zxfqcalving(i) + &
1146     fqcalving(i, nsrf)*pctsrf(i, nsrf)
1147     s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)
1148     s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)
1149     s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)
1150     s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)
1151     s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)
1152     s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)
1153     s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)
1154     s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)
1155     s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)
1156     s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)
1157     ENDDO
1158     ENDDO
1159    
1160     ! Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1161    
1162     DO nsrf = 1, nbsrf
1163     DO i = 1, klon
1164 guez 7 IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
1165 guez 3
1166 guez 7 IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)
1167     IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)
1168     IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)
1169     IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)
1170     IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)
1171     IF (pctsrf(i, nsrf) < epsfra) &
1172 guez 3 fqcalving(i, nsrf) = zxfqcalving(i)
1173 guez 7 IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf)=s_pblh(i)
1174     IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf)=s_lcl(i)
1175     IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf)=s_capCL(i)
1176     IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf)=s_oliqCL(i)
1177     IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf)=s_cteiCL(i)
1178     IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf)=s_pblT(i)
1179     IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf)=s_therm(i)
1180     IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf)=s_trmb1(i)
1181     IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf)=s_trmb2(i)
1182     IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf)=s_trmb3(i)
1183 guez 3 ENDDO
1184     ENDDO
1185    
1186     ! Calculer la derive du flux infrarouge
1187    
1188     DO i = 1, klon
1189     dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3
1190     ENDDO
1191    
1192     ! Appeler la convection (au choix)
1193    
1194     DO k = 1, llm
1195     DO i = 1, klon
1196     conv_q(i, k) = d_q_dyn(i, k) &
1197 guez 12 + d_q_vdf(i, k)/pdtphys
1198 guez 3 conv_t(i, k) = d_t_dyn(i, k) &
1199 guez 12 + d_t_vdf(i, k)/pdtphys
1200 guez 3 ENDDO
1201     ENDDO
1202     IF (check) THEN
1203     za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1204 guez 12 print *, "avantcon=", za
1205 guez 3 ENDIF
1206     zx_ajustq = .FALSE.
1207     IF (iflag_con == 2) zx_ajustq=.TRUE.
1208     IF (zx_ajustq) THEN
1209     DO i = 1, klon
1210     z_avant(i) = 0.0
1211     ENDDO
1212     DO k = 1, llm
1213     DO i = 1, klon
1214     z_avant(i) = z_avant(i) + (q_seri(i, k)+ql_seri(i, k)) &
1215 guez 17 *zmasse(i, k)
1216 guez 3 ENDDO
1217     ENDDO
1218     ENDIF
1219     IF (iflag_con == 1) THEN
1220     stop 'reactiver le call conlmd dans physiq.F'
1221     ELSE IF (iflag_con == 2) THEN
1222 guez 12 CALL conflx(pdtphys, paprs, pplay, t_seri, q_seri, &
1223 guez 3 conv_t, conv_q, zxfluxq(1, 1), omega, &
1224     d_t_con, d_q_con, rain_con, snow_con, &
1225     pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
1226     kcbot, kctop, kdtop, pmflxr, pmflxs)
1227     WHERE (rain_con < 0.) rain_con = 0.
1228     WHERE (snow_con < 0.) snow_con = 0.
1229     DO i = 1, klon
1230     ibas_con(i) = llm+1 - kcbot(i)
1231     itop_con(i) = llm+1 - kctop(i)
1232     ENDDO
1233     ELSE IF (iflag_con >= 3) THEN
1234     ! nb of tracers for the KE convection:
1235     ! MAF la partie traceurs est faite dans phytrac
1236     ! on met ntra=1 pour limiter les appels mais on peut
1237     ! supprimer les calculs / ftra.
1238     ntra = 1
1239     ! Schema de convection modularise et vectorise:
1240     ! (driver commun aux versions 3 et 4)
1241    
1242     IF (ok_cvl) THEN ! new driver for convectL
1243 guez 18 CALL concvl(iflag_con, pdtphys, paprs, pplay, t_seri, q_seri, &
1244 guez 3 u_seri, v_seri, tr_seri, ntra, &
1245     ema_work1, ema_work2, &
1246     d_t_con, d_q_con, d_u_con, d_v_con, d_tr, &
1247     rain_con, snow_con, ibas_con, itop_con, &
1248     upwd, dnwd, dnwd0, &
1249     Ma, cape, tvp, iflagctrl, &
1250     pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, wd, &
1251     pmflxr, pmflxs, &
1252     da, phi, mp)
1253    
1254     clwcon0=qcondc
1255 guez 13 pmfu=upwd+dnwd
1256 guez 3 ELSE ! ok_cvl
1257     ! MAF conema3 ne contient pas les traceurs
1258 guez 17 CALL conema3 (pdtphys, paprs, pplay, t_seri, q_seri, &
1259 guez 3 u_seri, v_seri, tr_seri, ntra, &
1260     ema_work1, ema_work2, &
1261     d_t_con, d_q_con, d_u_con, d_v_con, d_tr, &
1262     rain_con, snow_con, ibas_con, itop_con, &
1263     upwd, dnwd, dnwd0, bas, top, &
1264     Ma, cape, tvp, rflag, &
1265     pbase &
1266     , bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr &
1267     , clwcon0)
1268     ENDIF ! ok_cvl
1269    
1270     IF (.NOT. ok_gust) THEN
1271     do i = 1, klon
1272     wd(i)=0.0
1273     enddo
1274     ENDIF
1275    
1276     ! Calcul des proprietes des nuages convectifs
1277    
1278     DO k = 1, llm
1279     DO i = 1, klon
1280     zx_t = t_seri(i, k)
1281     IF (thermcep) THEN
1282     zdelta = MAX(0., SIGN(1., rtt-zx_t))
1283     zx_qs = r2es * FOEEW(zx_t, zdelta)/pplay(i, k)
1284     zx_qs = MIN(0.5, zx_qs)
1285     zcor = 1./(1.-retv*zx_qs)
1286     zx_qs = zx_qs*zcor
1287     ELSE
1288 guez 7 IF (zx_t < t_coup) THEN
1289 guez 3 zx_qs = qsats(zx_t)/pplay(i, k)
1290     ELSE
1291     zx_qs = qsatl(zx_t)/pplay(i, k)
1292     ENDIF
1293     ENDIF
1294     zqsat(i, k)=zx_qs
1295     ENDDO
1296     ENDDO
1297    
1298     ! calcul des proprietes des nuages convectifs
1299 guez 13 clwcon0=fact_cldcon*clwcon0
1300 guez 3 call clouds_gno &
1301     (klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, rnebcon0)
1302     ELSE
1303 guez 12 print *, "iflag_con non-prevu", iflag_con
1304 guez 3 stop 1
1305     ENDIF
1306    
1307     DO k = 1, llm
1308     DO i = 1, klon
1309     t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
1310     q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
1311     u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
1312     v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
1313     ENDDO
1314     ENDDO
1315    
1316     IF (if_ebil >= 2) THEN
1317     ztit='after convect'
1318 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1319 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1320 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1321     call diagphy(airephy, ztit, ip_ebil &
1322     , zero_v, zero_v, zero_v, zero_v, zero_v &
1323     , zero_v, rain_con, snow_con, ztsol &
1324     , d_h_vcol, d_qt, d_ec &
1325     , fs_bound, fq_bound )
1326     END IF
1327    
1328     IF (check) THEN
1329     za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1330 guez 12 print *,"aprescon=", za
1331 guez 3 zx_t = 0.0
1332     za = 0.0
1333     DO i = 1, klon
1334     za = za + airephy(i)/REAL(klon)
1335     zx_t = zx_t + (rain_con(i)+ &
1336     snow_con(i))*airephy(i)/REAL(klon)
1337     ENDDO
1338 guez 12 zx_t = zx_t/za*pdtphys
1339     print *,"Precip=", zx_t
1340 guez 3 ENDIF
1341     IF (zx_ajustq) THEN
1342     DO i = 1, klon
1343     z_apres(i) = 0.0
1344     ENDDO
1345     DO k = 1, llm
1346     DO i = 1, klon
1347     z_apres(i) = z_apres(i) + (q_seri(i, k)+ql_seri(i, k)) &
1348 guez 17 *zmasse(i, k)
1349 guez 3 ENDDO
1350     ENDDO
1351     DO i = 1, klon
1352 guez 12 z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*pdtphys) &
1353 guez 3 /z_apres(i)
1354     ENDDO
1355     DO k = 1, llm
1356     DO i = 1, klon
1357     IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
1358 guez 7 z_factor(i) < (1.0-1.0E-08)) THEN
1359 guez 3 q_seri(i, k) = q_seri(i, k) * z_factor(i)
1360     ENDIF
1361     ENDDO
1362     ENDDO
1363     ENDIF
1364     zx_ajustq=.FALSE.
1365    
1366     ! Convection seche (thermiques ou ajustement)
1367    
1368 guez 13 d_t_ajs=0.
1369     d_u_ajs=0.
1370     d_v_ajs=0.
1371     d_q_ajs=0.
1372     fm_therm=0.
1373     entr_therm=0.
1374 guez 3
1375 guez 12 IF(prt_level>9)print *, &
1376 guez 3 'AVANT LA CONVECTION SECHE, iflag_thermals=' &
1377     , iflag_thermals, ' nsplit_thermals=', nsplit_thermals
1378 guez 7 if(iflag_thermals < 0) then
1379 guez 3 ! Rien
1380 guez 12 IF(prt_level>9)print *,'pas de convection'
1381 guez 3 else if(iflag_thermals == 0) then
1382     ! Ajustement sec
1383 guez 12 IF(prt_level>9)print *,'ajsec'
1384 guez 3 CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
1385 guez 13 t_seri = t_seri + d_t_ajs
1386     q_seri = q_seri + d_q_ajs
1387 guez 3 else
1388     ! Thermiques
1389 guez 12 IF(prt_level>9)print *,'JUSTE AVANT, iflag_thermals=' &
1390 guez 3 , iflag_thermals, ' nsplit_thermals=', nsplit_thermals
1391     call calltherm(pdtphys &
1392     , pplay, paprs, pphi &
1393     , u_seri, v_seri, t_seri, q_seri &
1394     , d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs &
1395     , fm_therm, entr_therm)
1396     endif
1397    
1398     IF (if_ebil >= 2) THEN
1399     ztit='after dry_adjust'
1400 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1401 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1402 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1403     END IF
1404    
1405     ! Caclul des ratqs
1406    
1407     ! ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
1408     ! on ecrase le tableau ratqsc calcule par clouds_gno
1409     if (iflag_cldcon == 1) then
1410     do k=1, llm
1411     do i=1, klon
1412     if(ptconv(i, k)) then
1413     ratqsc(i, k)=ratqsbas &
1414     +fact_cldcon*(q_seri(i, 1)-q_seri(i, k))/q_seri(i, k)
1415     else
1416     ratqsc(i, k)=0.
1417     endif
1418     enddo
1419     enddo
1420     endif
1421    
1422     ! ratqs stables
1423     do k=1, llm
1424     do i=1, klon
1425     ratqss(i, k)=ratqsbas+(ratqshaut-ratqsbas)* &
1426     min((paprs(i, 1)-pplay(i, k))/(paprs(i, 1)-30000.), 1.)
1427     enddo
1428     enddo
1429    
1430     ! ratqs final
1431     if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then
1432     ! les ratqs sont une conbinaison de ratqss et ratqsc
1433     ! ratqs final
1434     ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1435     ! relaxation des ratqs
1436     facteur=exp(-pdtphys*facttemps)
1437 guez 13 ratqs=max(ratqs*facteur, ratqss)
1438     ratqs=max(ratqs, ratqsc)
1439 guez 3 else
1440     ! on ne prend que le ratqs stable pour fisrtilp
1441 guez 13 ratqs=ratqss
1442 guez 3 endif
1443    
1444     ! Appeler le processus de condensation a grande echelle
1445     ! et le processus de precipitation
1446 guez 12 CALL fisrtilp(pdtphys, paprs, pplay, &
1447 guez 3 t_seri, q_seri, ptconv, ratqs, &
1448     d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
1449     rain_lsc, snow_lsc, &
1450     pfrac_impa, pfrac_nucl, pfrac_1nucl, &
1451     frac_impa, frac_nucl, &
1452     prfl, psfl, rhcl)
1453    
1454     WHERE (rain_lsc < 0) rain_lsc = 0.
1455     WHERE (snow_lsc < 0) snow_lsc = 0.
1456     DO k = 1, llm
1457     DO i = 1, klon
1458     t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
1459     q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
1460     ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
1461     cldfra(i, k) = rneb(i, k)
1462     IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
1463     ENDDO
1464     ENDDO
1465     IF (check) THEN
1466     za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1467 guez 12 print *,"apresilp=", za
1468 guez 3 zx_t = 0.0
1469     za = 0.0
1470     DO i = 1, klon
1471     za = za + airephy(i)/REAL(klon)
1472     zx_t = zx_t + (rain_lsc(i) &
1473     + snow_lsc(i))*airephy(i)/REAL(klon)
1474     ENDDO
1475 guez 12 zx_t = zx_t/za*pdtphys
1476     print *,"Precip=", zx_t
1477 guez 3 ENDIF
1478    
1479     IF (if_ebil >= 2) THEN
1480     ztit='after fisrt'
1481 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1482 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1483 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1484     call diagphy(airephy, ztit, ip_ebil &
1485     , zero_v, zero_v, zero_v, zero_v, zero_v &
1486     , zero_v, rain_lsc, snow_lsc, ztsol &
1487     , d_h_vcol, d_qt, d_ec &
1488     , fs_bound, fq_bound )
1489     END IF
1490    
1491     ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1492    
1493     ! 1. NUAGES CONVECTIFS
1494    
1495     IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
1496     snow_tiedtke=0.
1497     if (iflag_cldcon == -1) then
1498     rain_tiedtke=rain_con
1499     else
1500     rain_tiedtke=0.
1501     do k=1, llm
1502     do i=1, klon
1503 guez 7 if (d_q_con(i, k) < 0.) then
1504 guez 3 rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i, k)/pdtphys &
1505 guez 17 *zmasse(i, k)
1506 guez 3 endif
1507     enddo
1508     enddo
1509     endif
1510    
1511     ! Nuages diagnostiques pour Tiedtke
1512     CALL diagcld1(paprs, pplay, &
1513     rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &
1514     diafra, dialiq)
1515     DO k = 1, llm
1516     DO i = 1, klon
1517     IF (diafra(i, k).GT.cldfra(i, k)) THEN
1518     cldliq(i, k) = dialiq(i, k)
1519     cldfra(i, k) = diafra(i, k)
1520     ENDIF
1521     ENDDO
1522     ENDDO
1523    
1524     ELSE IF (iflag_cldcon == 3) THEN
1525 guez 7 ! On prend pour les nuages convectifs le max du calcul de la
1526     ! convection et du calcul du pas de temps précédent diminué d'un facteur
1527     ! facttemps
1528 guez 3 facteur = pdtphys *facttemps
1529     do k=1, llm
1530     do i=1, klon
1531     rnebcon(i, k)=rnebcon(i, k)*facteur
1532     if (rnebcon0(i, k)*clwcon0(i, k).gt.rnebcon(i, k)*clwcon(i, k)) &
1533     then
1534     rnebcon(i, k)=rnebcon0(i, k)
1535     clwcon(i, k)=clwcon0(i, k)
1536     endif
1537     enddo
1538     enddo
1539    
1540     ! On prend la somme des fractions nuageuses et des contenus en eau
1541 guez 13 cldfra=min(max(cldfra, rnebcon), 1.)
1542     cldliq=cldliq+rnebcon*clwcon
1543 guez 3
1544     ENDIF
1545    
1546     ! 2. NUAGES STARTIFORMES
1547    
1548     IF (ok_stratus) THEN
1549     CALL diagcld2(paprs, pplay, t_seri, q_seri, diafra, dialiq)
1550     DO k = 1, llm
1551     DO i = 1, klon
1552     IF (diafra(i, k).GT.cldfra(i, k)) THEN
1553     cldliq(i, k) = dialiq(i, k)
1554     cldfra(i, k) = diafra(i, k)
1555     ENDIF
1556     ENDDO
1557     ENDDO
1558     ENDIF
1559    
1560     ! Precipitation totale
1561    
1562     DO i = 1, klon
1563     rain_fall(i) = rain_con(i) + rain_lsc(i)
1564     snow_fall(i) = snow_con(i) + snow_lsc(i)
1565     ENDDO
1566    
1567     IF (if_ebil >= 2) THEN
1568     ztit="after diagcld"
1569 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1570 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1571 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1572     END IF
1573    
1574     ! Calculer l'humidite relative pour diagnostique
1575    
1576     DO k = 1, llm
1577     DO i = 1, klon
1578     zx_t = t_seri(i, k)
1579     IF (thermcep) THEN
1580     zdelta = MAX(0., SIGN(1., rtt-zx_t))
1581     zx_qs = r2es * FOEEW(zx_t, zdelta)/pplay(i, k)
1582     zx_qs = MIN(0.5, zx_qs)
1583     zcor = 1./(1.-retv*zx_qs)
1584     zx_qs = zx_qs*zcor
1585     ELSE
1586 guez 7 IF (zx_t < t_coup) THEN
1587 guez 3 zx_qs = qsats(zx_t)/pplay(i, k)
1588     ELSE
1589     zx_qs = qsatl(zx_t)/pplay(i, k)
1590     ENDIF
1591     ENDIF
1592     zx_rh(i, k) = q_seri(i, k)/zx_qs
1593     zqsat(i, k)=zx_qs
1594     ENDDO
1595     ENDDO
1596     !jq - introduce the aerosol direct and first indirect radiative forcings
1597     !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
1598     IF (ok_ade.OR.ok_aie) THEN
1599     ! Get sulfate aerosol distribution
1600 guez 7 CALL readsulfate(rdayvrai, firstcal, sulfate)
1601     CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)
1602 guez 3
1603     ! Calculate aerosol optical properties (Olivier Boucher)
1604     CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, &
1605     tau_ae, piz_ae, cg_ae, aerindex)
1606     ELSE
1607     tau_ae(:, :, :)=0.0
1608     piz_ae(:, :, :)=0.0
1609     cg_ae(:, :, :)=0.0
1610     ENDIF
1611    
1612     ! Calculer les parametres optiques des nuages et quelques
1613     ! parametres pour diagnostiques:
1614    
1615     if (ok_newmicro) then
1616     CALL newmicro (paprs, pplay, ok_newmicro, &
1617     t_seri, cldliq, cldfra, cldtau, cldemi, &
1618     cldh, cldl, cldm, cldt, cldq, &
1619     flwp, fiwp, flwc, fiwc, &
1620     ok_aie, &
1621     sulfate, sulfate_pi, &
1622     bl95_b0, bl95_b1, &
1623     cldtaupi, re, fl)
1624     else
1625     CALL nuage (paprs, pplay, &
1626     t_seri, cldliq, cldfra, cldtau, cldemi, &
1627     cldh, cldl, cldm, cldt, cldq, &
1628     ok_aie, &
1629     sulfate, sulfate_pi, &
1630     bl95_b0, bl95_b1, &
1631     cldtaupi, re, fl)
1632    
1633     endif
1634    
1635     ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1636    
1637     IF (MOD(itaprad, radpas) == 0) THEN
1638     DO i = 1, klon
1639     albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &
1640     + falbe(i, is_lic) * pctsrf(i, is_lic) &
1641     + falbe(i, is_ter) * pctsrf(i, is_ter) &
1642     + falbe(i, is_sic) * pctsrf(i, is_sic)
1643     albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &
1644     + falblw(i, is_lic) * pctsrf(i, is_lic) &
1645     + falblw(i, is_ter) * pctsrf(i, is_ter) &
1646     + falblw(i, is_sic) * pctsrf(i, is_sic)
1647     ENDDO
1648     ! nouveau rayonnement (compatible Arpege-IFS):
1649     CALL radlwsw(dist, rmu0, fract, &
1650     paprs, pplay, zxtsol, albsol, albsollw, t_seri, q_seri, &
1651     wo, &
1652     cldfra, cldemi, cldtau, &
1653     heat, heat0, cool, cool0, radsol, albpla, &
1654     topsw, toplw, solsw, sollw, &
1655     sollwdown, &
1656     topsw0, toplw0, solsw0, sollw0, &
1657     lwdn0, lwdn, lwup0, lwup, &
1658     swdn0, swdn, swup0, swup, &
1659     ok_ade, ok_aie, & ! new for aerosol radiative effects
1660     tau_ae, piz_ae, cg_ae, &
1661     topswad, solswad, &
1662     cldtaupi, &
1663     topswai, solswai)
1664     itaprad = 0
1665     ENDIF
1666     itaprad = itaprad + 1
1667    
1668     ! Ajouter la tendance des rayonnements (tous les pas)
1669    
1670     DO k = 1, llm
1671     DO i = 1, klon
1672     t_seri(i, k) = t_seri(i, k) &
1673 guez 12 + (heat(i, k)-cool(i, k)) * pdtphys/86400.
1674 guez 3 ENDDO
1675     ENDDO
1676    
1677     IF (if_ebil >= 2) THEN
1678     ztit='after rad'
1679 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1680 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1681 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1682     call diagphy(airephy, ztit, ip_ebil &
1683     , topsw, toplw, solsw, sollw, zero_v &
1684     , zero_v, zero_v, zero_v, ztsol &
1685     , d_h_vcol, d_qt, d_ec &
1686     , fs_bound, fq_bound )
1687     END IF
1688    
1689     ! Calculer l'hydrologie de la surface
1690    
1691     DO i = 1, klon
1692     zxqsurf(i) = 0.0
1693     zxsnow(i) = 0.0
1694     ENDDO
1695     DO nsrf = 1, nbsrf
1696     DO i = 1, klon
1697     zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)
1698     zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)
1699     ENDDO
1700     ENDDO
1701    
1702     ! Calculer le bilan du sol et la derive de temperature (couplage)
1703    
1704     DO i = 1, klon
1705     bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1706     ENDDO
1707    
1708 guez 13 !mod deb lott(jan95)
1709 guez 3 ! Appeler le programme de parametrisation de l'orographie
1710     ! a l'echelle sous-maille:
1711    
1712     IF (ok_orodr) THEN
1713     ! selection des points pour lesquels le shema est actif:
1714     igwd=0
1715     DO i=1, klon
1716     itest(i)=0
1717     IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1718     itest(i)=1
1719     igwd=igwd+1
1720     idx(igwd)=i
1721     ENDIF
1722     ENDDO
1723    
1724 guez 12 CALL drag_noro(klon, llm, pdtphys, paprs, pplay, &
1725 guez 3 zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1726     igwd, idx, itest, &
1727     t_seri, u_seri, v_seri, &
1728     zulow, zvlow, zustrdr, zvstrdr, &
1729     d_t_oro, d_u_oro, d_v_oro)
1730    
1731     ! ajout des tendances
1732     DO k = 1, llm
1733     DO i = 1, klon
1734     t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
1735     u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
1736     v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
1737     ENDDO
1738     ENDDO
1739 guez 13 ENDIF
1740 guez 3
1741     IF (ok_orolf) THEN
1742    
1743     ! selection des points pour lesquels le shema est actif:
1744     igwd=0
1745     DO i=1, klon
1746     itest(i)=0
1747     IF ((zpic(i)-zmea(i)).GT.100.) THEN
1748     itest(i)=1
1749     igwd=igwd+1
1750     idx(igwd)=i
1751     ENDIF
1752     ENDDO
1753    
1754 guez 12 CALL lift_noro(klon, llm, pdtphys, paprs, pplay, &
1755 guez 3 rlat, zmea, zstd, zpic, &
1756     itest, &
1757     t_seri, u_seri, v_seri, &
1758     zulow, zvlow, zustrli, zvstrli, &
1759     d_t_lif, d_u_lif, d_v_lif)
1760    
1761     ! ajout des tendances
1762     DO k = 1, llm
1763     DO i = 1, klon
1764     t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
1765     u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
1766     v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
1767     ENDDO
1768     ENDDO
1769    
1770     ENDIF ! fin de test sur ok_orolf
1771    
1772     ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
1773    
1774     DO i = 1, klon
1775     zustrph(i)=0.
1776     zvstrph(i)=0.
1777     ENDDO
1778     DO k = 1, llm
1779     DO i = 1, klon
1780 guez 17 zustrph(i)=zustrph(i)+(u_seri(i, k)-u(i, k))/pdtphys* zmasse(i, k)
1781     zvstrph(i)=zvstrph(i)+(v_seri(i, k)-v(i, k))/pdtphys* zmasse(i, k)
1782 guez 3 ENDDO
1783     ENDDO
1784    
1785     !IM calcul composantes axiales du moment angulaire et couple des montagnes
1786    
1787 guez 17 CALL aaam_bud(27, klon, llm, gmtime, &
1788 guez 3 ra, rg, romega, &
1789     rlat, rlon, pphis, &
1790     zustrdr, zustrli, zustrph, &
1791     zvstrdr, zvstrli, zvstrph, &
1792     paprs, u, v, &
1793     aam, torsfc)
1794    
1795     IF (if_ebil >= 2) THEN
1796     ztit='after orography'
1797 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1798 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1799 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1800     END IF
1801    
1802     !AA Installation de l'interface online-offline pour traceurs
1803    
1804     ! Calcul des tendances traceurs
1805    
1806 guez 7 call phytrac(rnpb, itap, lmt_pas, julien, gmtime, firstcal, lafin, nq-2, &
1807 guez 12 pdtphys, u, v, t, paprs, pplay, pmfu, pmfd, pen_u, pde_u, pen_d, &
1808     pde_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, &
1809 guez 20 frac_impa, frac_nucl, pphis, pphi, albsol, rhcl, cldfra, &
1810 guez 12 rneb, diafra, cldliq, itop_con, ibas_con, pmflxr, pmflxs, prfl, &
1811 guez 17 psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)
1812 guez 3
1813     IF (offline) THEN
1814    
1815     print*, 'Attention on met a 0 les thermiques pour phystoke'
1816     call phystokenc(pdtphys, rlon, rlat, &
1817     t, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
1818     fm_therm, entr_therm, &
1819     ycoefh, yu1, yv1, ftsol, pctsrf, &
1820     frac_impa, frac_nucl, &
1821 guez 12 pphis, airephy, pdtphys, itap)
1822 guez 3
1823     ENDIF
1824    
1825     ! Calculer le transport de l'eau et de l'energie (diagnostique)
1826    
1827     CALL transp (paprs, zxtsol, &
1828     t_seri, q_seri, u_seri, v_seri, zphi, &
1829     ve, vq, ue, uq)
1830    
1831     !IM diag. bilKP
1832    
1833     CALL transp_lay (paprs, zxtsol, &
1834     t_seri, q_seri, u_seri, v_seri, zphi, &
1835     ve_lay, vq_lay, ue_lay, uq_lay)
1836    
1837     ! Accumuler les variables a stocker dans les fichiers histoire:
1838    
1839     !+jld ec_conser
1840     DO k = 1, llm
1841     DO i = 1, klon
1842     ZRCPD = RCPD*(1.0+RVTMP2*q_seri(i, k))
1843     d_t_ec(i, k)=0.5/ZRCPD &
1844     *(u(i, k)**2+v(i, k)**2-u_seri(i, k)**2-v_seri(i, k)**2)
1845     t_seri(i, k)=t_seri(i, k)+d_t_ec(i, k)
1846 guez 12 d_t_ec(i, k) = d_t_ec(i, k)/pdtphys
1847 guez 3 END DO
1848     END DO
1849     !-jld ec_conser
1850     IF (if_ebil >= 1) THEN
1851     ztit='after physic'
1852 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, pdtphys &
1853 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1854 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1855     ! Comme les tendances de la physique sont ajoute dans la dynamique,
1856     ! on devrait avoir que la variation d'entalpie par la dynamique
1857     ! est egale a la variation de la physique au pas de temps precedent.
1858     ! Donc la somme de ces 2 variations devrait etre nulle.
1859     call diagphy(airephy, ztit, ip_ebil &
1860     , topsw, toplw, solsw, sollw, sens &
1861     , evap, rain_fall, snow_fall, ztsol &
1862     , d_h_vcol, d_qt, d_ec &
1863     , fs_bound, fq_bound )
1864    
1865     d_h_vcol_phy=d_h_vcol
1866    
1867     END IF
1868    
1869     ! SORTIES
1870    
1871     !cc prw = eau precipitable
1872     DO i = 1, klon
1873     prw(i) = 0.
1874     DO k = 1, llm
1875 guez 17 prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)
1876 guez 3 ENDDO
1877     ENDDO
1878    
1879     ! Convertir les incrementations en tendances
1880    
1881     DO k = 1, llm
1882     DO i = 1, klon
1883 guez 12 d_u(i, k) = ( u_seri(i, k) - u(i, k) ) / pdtphys
1884     d_v(i, k) = ( v_seri(i, k) - v(i, k) ) / pdtphys
1885     d_t(i, k) = ( t_seri(i, k)-t(i, k) ) / pdtphys
1886     d_qx(i, k, ivap) = ( q_seri(i, k) - qx(i, k, ivap) ) / pdtphys
1887     d_qx(i, k, iliq) = ( ql_seri(i, k) - qx(i, k, iliq) ) / pdtphys
1888 guez 3 ENDDO
1889     ENDDO
1890    
1891     IF (nq >= 3) THEN
1892     DO iq = 3, nq
1893     DO k = 1, llm
1894     DO i = 1, klon
1895 guez 13 d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / pdtphys
1896 guez 3 ENDDO
1897     ENDDO
1898     ENDDO
1899     ENDIF
1900    
1901     ! Sauvegarder les valeurs de t et q a la fin de la physique:
1902     DO k = 1, llm
1903     DO i = 1, klon
1904     t_ancien(i, k) = t_seri(i, k)
1905     q_ancien(i, k) = q_seri(i, k)
1906     ENDDO
1907     ENDDO
1908    
1909     ! Ecriture des sorties
1910     call write_histhf
1911     call write_histday
1912     call write_histins
1913    
1914     ! Si c'est la fin, il faut conserver l'etat de redemarrage
1915     IF (lafin) THEN
1916     itau_phy = itau_phy + itap
1917 guez 13 CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, &
1918 guez 12 ftsoil, tslab, seaice, fqsurf, qsol, &
1919 guez 3 fsnow, falbe, falblw, fevap, rain_fall, snow_fall, &
1920     solsw, sollwdown, dlw, &
1921     radsol, frugs, agesno, &
1922 guez 13 zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1923 guez 3 t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
1924     ENDIF
1925    
1926     contains
1927    
1928 guez 15 subroutine write_histday
1929 guez 3
1930 guez 15 use grid_change, only: gr_phy_write_3d
1931 guez 17 integer itau_w ! pas de temps ecriture
1932 guez 3
1933 guez 15 !------------------------------------------------
1934 guez 3
1935     if (ok_journe) THEN
1936     itau_w = itau_phy + itap
1937 guez 17 if (nq <= 4) then
1938     call histwrite(nid_day, "Sigma_O3_Royer", itau_w, &
1939     gr_phy_write_3d(wo) * 1e3)
1940     ! (convert "wo" from kDU to DU)
1941     end if
1942 guez 3 if (ok_sync) then
1943     call histsync(nid_day)
1944     endif
1945     ENDIF
1946    
1947     End subroutine write_histday
1948    
1949     !****************************
1950    
1951     subroutine write_histhf
1952    
1953     ! From phylmd/write_histhf.h, v 1.5 2005/05/25 13:10:09
1954    
1955 guez 17 !------------------------------------------------
1956 guez 3
1957     call write_histhf3d
1958    
1959     IF (ok_sync) THEN
1960     call histsync(nid_hf)
1961     ENDIF
1962    
1963     end subroutine write_histhf
1964    
1965     !***************************************************************
1966    
1967     subroutine write_histins
1968    
1969     ! From phylmd/write_histins.h, v 1.2 2005/05/25 13:10:09
1970    
1971     real zout
1972 guez 17 integer itau_w ! pas de temps ecriture
1973 guez 3
1974     !--------------------------------------------------
1975    
1976     IF (ok_instan) THEN
1977     ! Champs 2D:
1978    
1979 guez 12 zsto = pdtphys * ecrit_ins
1980     zout = pdtphys * ecrit_ins
1981 guez 3 itau_w = itau_phy + itap
1982    
1983     i = NINT(zout/zsto)
1984     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), pphis, zx_tmp_2d)
1985 guez 15 CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)
1986 guez 3
1987     i = NINT(zout/zsto)
1988     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), airephy, zx_tmp_2d)
1989 guez 15 CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)
1990 guez 3
1991     DO i = 1, klon
1992     zx_tmp_fi2d(i) = paprs(i, 1)
1993     ENDDO
1994     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
1995 guez 15 CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d)
1996 guez 3
1997     DO i = 1, klon
1998     zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
1999     ENDDO
2000     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2001 guez 15 CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)
2002 guez 3
2003     DO i = 1, klon
2004     zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
2005     ENDDO
2006     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2007 guez 15 CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)
2008 guez 3
2009     DO i = 1, klon
2010     zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
2011     ENDDO
2012     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2013 guez 15 CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)
2014 guez 3
2015     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zxtsol, zx_tmp_2d)
2016 guez 15 CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)
2017 guez 3 !ccIM
2018     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zt2m, zx_tmp_2d)
2019 guez 15 CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)
2020 guez 3
2021     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zq2m, zx_tmp_2d)
2022 guez 15 CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)
2023 guez 3
2024     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zu10m, zx_tmp_2d)
2025 guez 15 CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)
2026 guez 3
2027     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zv10m, zx_tmp_2d)
2028 guez 15 CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)
2029 guez 3
2030     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), snow_fall, zx_tmp_2d)
2031 guez 15 CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)
2032 guez 3
2033     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), cdragm, zx_tmp_2d)
2034 guez 15 CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)
2035 guez 3
2036     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), cdragh, zx_tmp_2d)
2037 guez 15 CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)
2038 guez 3
2039     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), toplw, zx_tmp_2d)
2040 guez 15 CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)
2041 guez 3
2042     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), evap, zx_tmp_2d)
2043 guez 15 CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)
2044 guez 3
2045     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), solsw, zx_tmp_2d)
2046 guez 15 CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)
2047 guez 3
2048     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sollw, zx_tmp_2d)
2049 guez 15 CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)
2050 guez 3
2051     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sollwdown, zx_tmp_2d)
2052 guez 15 CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)
2053 guez 3
2054     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), bils, zx_tmp_2d)
2055 guez 15 CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)
2056 guez 3
2057     zx_tmp_fi2d(1:klon)=-1*sens(1:klon)
2058     ! CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sens, zx_tmp_2d)
2059     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2060 guez 15 CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)
2061 guez 3
2062     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), fder, zx_tmp_2d)
2063 guez 15 CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)
2064 guez 3
2065     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_oce), zx_tmp_2d)
2066 guez 15 CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)
2067 guez 3
2068     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_ter), zx_tmp_2d)
2069 guez 15 CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)
2070 guez 3
2071     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_lic), zx_tmp_2d)
2072 guez 15 CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)
2073 guez 3
2074     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_sic), zx_tmp_2d)
2075 guez 15 CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)
2076 guez 3
2077     DO nsrf = 1, nbsrf
2078     !XXX
2079     zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)*100.
2080     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2081     CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &
2082 guez 15 zx_tmp_2d)
2083 guez 3
2084     zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
2085     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2086     CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &
2087 guez 15 zx_tmp_2d)
2088 guez 3
2089     zx_tmp_fi2d(1 : klon) = fluxt( 1 : klon, 1, nsrf)
2090     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2091     CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &
2092 guez 15 zx_tmp_2d)
2093 guez 3
2094     zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, nsrf)
2095     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2096     CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &
2097 guez 15 zx_tmp_2d)
2098 guez 3
2099     zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)
2100     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2101     CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &
2102 guez 15 zx_tmp_2d)
2103 guez 3
2104     zx_tmp_fi2d(1 : klon) = fluxu( 1 : klon, 1, nsrf)
2105     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2106     CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &
2107 guez 15 zx_tmp_2d)
2108 guez 3
2109     zx_tmp_fi2d(1 : klon) = fluxv( 1 : klon, 1, nsrf)
2110     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2111     CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &
2112 guez 15 zx_tmp_2d)
2113 guez 3
2114     zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)
2115     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2116     CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &
2117 guez 15 zx_tmp_2d)
2118 guez 3
2119     zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
2120     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2121     CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &
2122 guez 15 zx_tmp_2d)
2123 guez 3
2124     END DO
2125     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), albsol, zx_tmp_2d)
2126 guez 15 CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)
2127 guez 3 CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), albsollw, zx_tmp_2d)
2128 guez 15 CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)
2129 guez 3
2130     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zxrugs, zx_tmp_2d)
2131 guez 15 CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)
2132 guez 3
2133     !IM cf. AM 081204 BEG
2134    
2135     !HBTM2
2136    
2137     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_pblh, zx_tmp_2d)
2138 guez 15 CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d)
2139 guez 3
2140     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_pblt, zx_tmp_2d)
2141 guez 15 CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d)
2142 guez 3
2143     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_lcl, zx_tmp_2d)
2144 guez 15 CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d)
2145 guez 3
2146     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_capCL, zx_tmp_2d)
2147 guez 15 CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d)
2148 guez 3
2149     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_oliqCL, zx_tmp_2d)
2150 guez 15 CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d)
2151 guez 3
2152     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_cteiCL, zx_tmp_2d)
2153 guez 15 CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d)
2154 guez 3
2155     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_therm, zx_tmp_2d)
2156 guez 15 CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d)
2157 guez 3
2158     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb1, zx_tmp_2d)
2159 guez 15 CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d)
2160 guez 3
2161     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb2, zx_tmp_2d)
2162 guez 15 CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d)
2163 guez 3
2164     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb3, zx_tmp_2d)
2165 guez 15 CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)
2166 guez 3
2167     !IM cf. AM 081204 END
2168    
2169     ! Champs 3D:
2170    
2171     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), t_seri, zx_tmp_3d)
2172 guez 15 CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d)
2173 guez 3
2174     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), u_seri, zx_tmp_3d)
2175 guez 15 CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d)
2176 guez 3
2177     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), v_seri, zx_tmp_3d)
2178 guez 15 CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d)
2179 guez 3
2180     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), zphi, zx_tmp_3d)
2181 guez 15 CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d)
2182 guez 3
2183     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), pplay, zx_tmp_3d)
2184 guez 15 CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d)
2185 guez 3
2186     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), d_t_vdf, zx_tmp_3d)
2187 guez 15 CALL histwrite(nid_ins, "dtvdf", itau_w, zx_tmp_3d)
2188 guez 3
2189     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), d_q_vdf, zx_tmp_3d)
2190 guez 15 CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
2191 guez 3
2192     if (ok_sync) then
2193     call histsync(nid_ins)
2194     endif
2195     ENDIF
2196    
2197     end subroutine write_histins
2198    
2199     !****************************************************
2200    
2201     subroutine write_histhf3d
2202    
2203     ! From phylmd/write_histhf3d.h, v 1.2 2005/05/25 13:10:09
2204    
2205 guez 17 integer itau_w ! pas de temps ecriture
2206 guez 3
2207 guez 17 !-------------------------------------------------------
2208    
2209 guez 3 itau_w = itau_phy + itap
2210    
2211     ! Champs 3D:
2212    
2213     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), t_seri, zx_tmp_3d)
2214 guez 15 CALL histwrite(nid_hf3d, "temp", itau_w, zx_tmp_3d)
2215 guez 3
2216     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), qx(1, 1, ivap), zx_tmp_3d)
2217 guez 15 CALL histwrite(nid_hf3d, "ovap", itau_w, zx_tmp_3d)
2218 guez 3
2219     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), u_seri, zx_tmp_3d)
2220 guez 15 CALL histwrite(nid_hf3d, "vitu", itau_w, zx_tmp_3d)
2221 guez 3
2222     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), v_seri, zx_tmp_3d)
2223 guez 15 CALL histwrite(nid_hf3d, "vitv", itau_w, zx_tmp_3d)
2224 guez 3
2225 guez 6 if (nbtr >= 3) then
2226     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), tr_seri(1, 1, 3), &
2227     zx_tmp_3d)
2228 guez 15 CALL histwrite(nid_hf3d, "O3", itau_w, zx_tmp_3d)
2229 guez 6 end if
2230 guez 3
2231     if (ok_sync) then
2232     call histsync(nid_hf3d)
2233     endif
2234    
2235     end subroutine write_histhf3d
2236    
2237     END SUBROUTINE physiq
2238    
2239     end module physiq_m

  ViewVC Help
Powered by ViewVC 1.1.21