/[lmdze]/trunk/libf/phylmd/physiq.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/physiq.f90

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21