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

Annotation of /trunk/phylmd/physiq.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 20 - (hide annotations)
Wed Oct 15 16:19:57 2008 UTC (15 years, 7 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 76646 byte(s) 
Deleted argument "presnivs" of "physiq", "ini_histhf", "ini_histhf3d",
"ini_histday", "ini_histins", "ini_histrac", "phytrac". Access it from
"comvert" instead.

Replaced calls to NetCDF Fortran 77 interface by calls to Fortran 90
interface or to NetCDF95.

Procedure "gr_phy_write_3d" now works with a variable of arbitrary
size in the second dimension.

Annotated use statements with "only" clause.

Replaced calls to NetCDF interface version 2 by calls to Fortran 90
interface in "guide.f90" and "read_reanalyse.f".

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