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

Annotation of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 7 - (hide annotations)
Mon Mar 31 12:24:17 2008 UTC (16 years, 1 month ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 93342 byte(s) 
This revision is not in working order. Pending some moving of files.

Important changes. In the program "etat0_lim": ozone coefficients from
Mobidic are regridded in time instead of pressure ; consequences in
"etat0". In the program "gcm", ozone coefficients from Mobidic are
read once per day only for the current day and regridded in pressure ;
consequences in "o3_chem_m", "regr_pr_coefoz", "phytrac" and
"regr_pr_comb_coefoz_m".

NetCDF95 is a library and does not export NetCDF.

New variables "nag_gl_options", "nag_fcalls_options" and
"nag_cross_options" in "nag_tools.mk".

"check_coefoz.jnl" rewritten entirely for new version of
"coefoz_LMDZ.nc".

Target "obj_etat0_lim" moved from "GNUmakefile" to "nag_rules.mk".

Added some "intent" attributes in "calfis", "clmain", "clqh",
"cltrac", "cltracrn", "cvltr", "ini_undefSTD", "moy_undefSTD",
"nflxtr", "phystokenc", "phytrac", "readsulfate", "readsulfate_preind"
and "undefSTD".

In "dynetat0", "dynredem0" and "gcm", "phis" has rank 2 instead of
1. "phis" has assumed shape in "dynredem0".

Added module containing "dynredem0". Changed some calls with NetCDF
Fortran 77 interface to calls with NetCDF95 interface.

Replaced calls to "ssum" by calls to "sum" in "inigeom".

In "make.sh", new option "-c" to change compiler.

In "aaam_bud", argument "rjour" deleted.

In "physiq": renamed some variables; deleted variable "xjour".

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