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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 17 - (hide annotations)
Tue Aug 5 13:31:32 2008 UTC (15 years, 9 months ago) by guez
File size: 76848 byte(s) 
Created rule for "compare_sampl_*" files in
"Documentation/Manuel_LMDZE.texfol/Graphiques/GNUmakefile".

Extracted "qcheck", "radiornpb", "minmaxqfi" into separate files.

Read pressure coordinate of ozone coefficients once per run instead of
every day.

Added some "intent" attributes.

Added argument "nq" to "ini_histday". Replaced calls to "gr_fi_ecrit"
by calls to "gr_phy_write_2d". "Sigma_O3_Royer" is written to
"histday.nc" only if "nq >= 4". Moved "ini_histrac" to module
"ini_hist".

Compute "zmasse" in "physiq", pass it to "phytrac".

Removed computations of "pftsol*" and "ppsrf*" in "phytrac".

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