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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 12 - (hide annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 9 months ago) by guez
File size: 92356 byte(s) 
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

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