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

Annotation of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 13 - (hide annotations)
Fri Jul 25 19:59:34 2008 UTC (15 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 91675 byte(s)
-- Minor change of behaviour:

"etat0" does not compute "rugsrel" nor "radpas". Deleted arguments
"radpas" and "rugsrel" of "phyredem". Deleted argument "rugsrel" of
"phyetat0". "startphy.nc" does not contain the variable "RUGSREL". In
"physiq", "rugoro" is set to 0 if not "ok_orodr". The whole program
"etat0_lim" does not use "clesphys2".

-- Minor modification of input/output:

Created subroutine "read_clesphys2". Variables of "clesphys2" are read
in "read_clesphys2" instead of "conf_gcm". "printflag" does not print
variables of "clesphys2".

-- Should not change any result at run time:

References to module "numer_rec" instead of individual modules of
"Numer_rec_Lionel".

Deleted argument "clesphy0" of "calfis", "physiq", "conf_gcm",
"leapfrog", "phyetat0". Deleted variable "clesphy0" in
"gcm". "phyetat0" does not modify variables of "clesphys2".

The program unit "gcm" does not modify "itau_phy".

Added some "intent" attributes.

"regr11_lint" does not call "polint".

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

  ViewVC Help
Powered by ViewVC 1.1.21