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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 18 - (hide annotations)
Thu Aug 7 12:29:13 2008 UTC (15 years, 9 months ago) by guez
File size: 76785 byte(s) 
In module "regr_pr", rewrote scanning of horizontal positions as a
single set of loops, using a mask.

Added some "intent" attributes.

In "dynredem0", replaced calls to Fortran 77 interface of NetCDF by
calls to NetCDF95. Removed calls to "nf_redef", regrouped all writing
operations. In "dynredem1", replaced some calls to Fortran 77
interface of NetCDF by calls to Fortran 90 interface.

Renamed variable "nqmax" to "nq_phys".

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