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

Annotation of /trunk/phylmd/physiq.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 16 - (hide annotations)
Fri Aug 1 15:37:00 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 78217 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.21