Context Navigation

physiq.F @ 1263

Last change on this file since 1263 was 1263, checked in by lguez, 15 years ago

1) Reactivated ability to read ozone (that was deactivated because of
dependency on version of IOIPSL). Added ability to read a pressure
coordinate in Pa in "regr_lat_time_climoz".

2) Added the ability to read a second ozone climatology, corresponding to
daylight ozone:

-- "read_climoz" is now an integer variable, instead of a logical
variable.

-- Added argument "read_climoz" to "phys_state_var_init",
"phys_output_open" and "regr_lat_time_climoz".

-- Created new variable "ozone_daylight" for "hist*.nc" output files.

-- Added a third dimension to variable "wo" in module
"phys_state_var_mod" and variable "POZON" in "radlwsw": index 1 for
average day-night ozone, index 2 for daylight ozone.

-- Added a fourth dimension to variables "o3_in", "o3_regr_lat" and
"o3_out" in "regr_lat_time_climoz": index 1 for average day-night
ozone, index 2 for daylight ozone.

-- In "physiq", moved call to "conf_phys" before call to
"phys_state_var_init". Thus, "conf_phys" is now inside the block "if
(first)" instead of "IF (debut)". There were definitions of "bl95_b0"
and "bl95_b1" that were useless because the variables were overwritten
by "conf_phys". Removed those definitions.

-- In "radlwsw", we pass the average day-night ozone to "LW_LMDAR4"
and the daylight ozone, if we have it, to "SW_LMDAR4" or
"SW_AEROAR4". If we do not have a specific field for daylight ozone
then "SW_LMDAR4" or "SW_AEROAR4" just get the average day-night ozone.

-- "regr_lat_time_climoz" now manages latitudes where the input ozone
field is missing at all levels (polar night).

-- Encapsulated "radlwsw" in a module.

3) Modifications to make sequential and parallel versions of
"create_etat0_limit" almost identical:

-- In "dyn3dpar/create_etat0_limit.F". No need to call
"phys_state_var_init", removed "use phys_state_var_mod" statement. No
need for "clesphys.h", removed "include" statement.

-- In "dyn3dpar/etat0_netcdf.F". Added argument "tau_ratqs" in call to
"conf_phys" (this bug was already corrected in "dyn3d"). Moved call to
"inifilr" after call to "infotrac_init" (as in "dyn3d").

4) Other peripheral modifications:

-- Added procedures "nf95_get_att" and "nf95_def_var_scalar" in
NetCDF95 interface. Overloaded "nf95_put_var" with three more
procedures: "nf95_put_var_FourByteReal", "nf95_put_var_FourByteInt",
"nf95_put_var_1D_FourByteInt".

-- Overloaded "regr1_step_av" with one more procedure:
"regr14_step_av". Overloaded "regr3_lint" with one more procedure:
"regr34_lint".

-- Corrected call to "Init_Phys_lmdz" in "dyn3d/create_etat0_limit.F":
the last argument should be an array, not a scalar.

-- Encapsulated "conf_phys" in a module.

-- Splitted module "regr_pr" into "regr_pr_av_m" and "regr_pr_int_m".

5) Tests:

This revision was compared to revision 1259, with optimization options
"debug" and "dev", parallelization options "none", "mpi", "omp" and
"mpi_omp", 1 and 2 MPI processes, 1 and 2 OpenMP threads, with the
compiler "FORTRAN90/SX Version 2.0 for SX-8". Both programs
"create_etat0_limit" and "gcm" were tested. In all cases,
parallelization does not change the results. With "read_climoz = 0" in
the ".def" files, the results of revision 1259 and of this revision
are the same.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 117.4 KB

Line
1	! $Id: physiq.F 1263 2009-11-17 13:00:14Z lguez $
2	!
3	c#define IO_DEBUG
4
5	SUBROUTINE physiq (nlon,nlev,
6	. debut,lafin,jD_cur, jH_cur,pdtphys,
7	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
8	. u,v,t,qx,
9	. flxmass_w,
10	. d_u, d_v, d_t, d_qx, d_ps
11	. , dudyn
12	. , PVteta)
13
14	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
15	$ histwrite, ju2ymds, ymds2ju, ioget_year_len
16	USE comgeomphy
17	USE phys_cal_mod
18	USE write_field_phy
19	USE dimphy
20	USE infotrac
21	USE mod_grid_phy_lmdz
22	USE mod_phys_lmdz_para
23	USE iophy
24	USE misc_mod, mydebug=>debug
25	USE vampir
26	USE pbl_surface_mod, ONLY : pbl_surface
27	USE change_srf_frac_mod
28	USE surface_data, ONLY : type_ocean, ok_veget
29	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
30	USE phys_state_var_mod ! Variables sauvegardees de la physique
31	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
32	USE phys_output_mod
33	use open_climoz_m, only: open_climoz ! ozone climatology from a file
34	use regr_pr_av_m, only: regr_pr_av
35	use netcdf95, only: nf95_close
36	use mod_phys_lmdz_mpi_data, only: is_mpi_root
37	USE aero_mod
38	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
39	use conf_phys_m, only: conf_phys
40	use radlwsw_m, only: radlwsw
41
42	IMPLICIT none
43	c======================================================================
44	c
45	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
46	c
47	c Objet: Moniteur general de la physique du modele
48	cAA Modifications quant aux traceurs :
49	cAA - uniformisation des parametrisations ds phytrac
50	cAA - stockage des moyennes des champs necessaires
51	cAA en mode traceur off-line
52	c======================================================================
53	c CLEFS CPP POUR LES IO
54	c =====================
55	c#define histmthNMC
56	c#define histISCCP
57	c======================================================================
58	c modif ( P. Le Van , 12/10/98 )
59	c
60	c Arguments:
61	c
62	c nlon----input-I-nombre de points horizontaux
63	c nlev----input-I-nombre de couches verticales, doit etre egale a klev
64	c debut---input-L-variable logique indiquant le premier passage
65	c lafin---input-L-variable logique indiquant le dernier passage
66	c jD_cur -R-jour courant a l'appel de la physique (jour julien)
67	c jH_cur -R-heure courante a l'appel de la physique (jour julien)
68	c pdtphys-input-R-pas d'integration pour la physique (seconde)
69	c paprs---input-R-pression pour chaque inter-couche (en Pa)
70	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
71	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
72	c pphis---input-R-geopotentiel du sol
73	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
74	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
75	c v-------input-R-vitesse Y (de S a N) en m/s
76	c t-------input-R-temperature (K)
77	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
78	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
79	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
80	c flxmass_w -input-R- flux de masse verticale
81	c d_u-----output-R-tendance physique de "u" (m/s/s)
82	c d_v-----output-R-tendance physique de "v" (m/s/s)
83	c d_t-----output-R-tendance physique de "t" (K/s)
84	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
85	c d_ps----output-R-tendance physique de la pression au sol
86	cIM
87	c PVteta--output-R-vorticite potentielle a des thetas constantes
88	c======================================================================
89	#include "dimensions.h"
90	integer jjmp1
91	parameter (jjmp1=jjm+1-1/jjm)
92	integer iip1
93	parameter (iip1=iim+1)
94
95	#include "regdim.h"
96	#include "indicesol.h"
97	#include "dimsoil.h"
98	#include "clesphys.h"
99	#include "control.h"
100	#include "temps.h"
101	#include "iniprint.h"
102	#include "thermcell.h"
103	c======================================================================
104	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
105	PARAMETER (ok_cvl=.TRUE.)
106	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
107	PARAMETER (ok_gust=.FALSE.)
108	integer iflag_radia ! active ou non le rayonnement (MPL)
109	save iflag_radia
110	c$OMP THREADPRIVATE(iflag_radia)
111	c======================================================================
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	c======================================================================
117	REAL amn, amx
118	INTEGER igout
119	c======================================================================
120	c Clef controlant l'activation du cycle diurne:
121	ccc LOGICAL cycle_diurne
122	ccc PARAMETER (cycle_diurne=.FALSE.)
123	c======================================================================
124	c Modele thermique du sol, a activer pour le cycle diurne:
125	ccc LOGICAL soil_model
126	ccc PARAMETER (soil_model=.FALSE.)
127	c======================================================================
128	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
129	c le calcul du rayonnement est celle apres la precipitation des nuages.
130	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
131	c la condensation et la precipitation. Cette cle augmente les impacts
132	c radiatifs des nuages.
133	ccc LOGICAL new_oliq
134	ccc PARAMETER (new_oliq=.FALSE.)
135	c======================================================================
136	c Clefs controlant deux parametrisations de l'orographie:
137	cc LOGICAL ok_orodr
138	ccc PARAMETER (ok_orodr=.FALSE.)
139	ccc LOGICAL ok_orolf
140	ccc PARAMETER (ok_orolf=.FALSE.)
141	c======================================================================
142	LOGICAL ok_journe ! sortir le fichier journalier
143	save ok_journe
144	c$OMP THREADPRIVATE(ok_journe)
145	c
146	LOGICAL ok_mensuel ! sortir le fichier mensuel
147	save ok_mensuel
148	c$OMP THREADPRIVATE(ok_mensuel)
149	c
150	LOGICAL ok_instan ! sortir le fichier instantane
151	save ok_instan
152	c$OMP THREADPRIVATE(ok_instan)
153	c
154	LOGICAL ok_LES ! sortir le fichier LES
155	save ok_LES
156	c$OMP THREADPRIVATE(ok_LES)
157	c
158	LOGICAL ok_region ! sortir le fichier regional
159	PARAMETER (ok_region=.FALSE.)
160	c======================================================================
161	real weak_inversion(klon),dthmin(klon)
162	real seuil_inversion
163	save seuil_inversion
164	c$OMP THREADPRIVATE(seuil_inversion)
165	integer iflag_ratqs
166	save iflag_ratqs
167	c$OMP THREADPRIVATE(iflag_ratqs)
168	REAL lambda_th(klon,klev),zz,znum,zden
169	REAL wmax_th(klon)
170	REAL zmax_th(klon)
171	REAL tau_overturning_th(klon)
172
173	integer lmax_th(klon)
174	integer limbas(klon)
175	real ratqscth(klon,klev)
176	real ratqsdiff(klon,klev)
177	real zqsatth(klon,klev)
178
179	c======================================================================
180	c
181	INTEGER ivap ! indice de traceurs pour vapeur d'eau
182	PARAMETER (ivap=1)
183	INTEGER iliq ! indice de traceurs pour eau liquide
184	PARAMETER (iliq=2)
185
186	c
187	c
188	c Variables argument:
189	c
190	INTEGER nlon
191	INTEGER nlev
192	REAL, intent(in):: jD_cur, jH_cur
193
194	REAL pdtphys
195	LOGICAL debut, lafin
196	REAL paprs(klon,klev+1)
197	REAL pplay(klon,klev)
198	REAL pphi(klon,klev)
199	REAL pphis(klon)
200	REAL presnivs(klev)
201	REAL znivsig(klev)
202	real pir
203
204	REAL u(klon,klev)
205	REAL v(klon,klev)
206	REAL t(klon,klev),theta(klon,klev)
207	REAL qx(klon,klev,nqtot)
208	REAL flxmass_w(klon,klev)
209	REAL omega(klon,klev) ! vitesse verticale en Pa/s
210	REAL d_u(klon,klev)
211	REAL d_v(klon,klev)
212	REAL d_t(klon,klev)
213	REAL d_qx(klon,klev,nqtot)
214	REAL d_ps(klon)
215	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
216	c
217	cIM Amip2 PV a theta constante
218	c
219	INTEGER nbteta
220	PARAMETER(nbteta=3)
221	CHARACTER*3 ctetaSTD(nbteta)
222	DATA ctetaSTD/'350','380','405'/
223	SAVE ctetaSTD
224	c$OMP THREADPRIVATE(ctetaSTD)
225	REAL rtetaSTD(nbteta)
226	DATA rtetaSTD/350., 380., 405./
227	SAVE rtetaSTD
228	c$OMP THREADPRIVATE(rtetaSTD)
229	c
230	REAL PVteta(klon,nbteta)
231	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
232	c
233	cMI Amip2 PV a theta constante
234
235	cym INTEGER klevp1, klevm1
236	cym PARAMETER(klevp1=klev+1,klevm1=klev-1)
237	cym#include "raddim.h"
238	c
239	c
240	cIM Amip2
241	c variables a une pression donnee
242	c
243	real rlevSTD(nlevSTD)
244	DATA rlevSTD/100000., 92500., 85000., 70000.,
245	.60000., 50000., 40000., 30000., 25000., 20000.,
246	.15000., 10000., 7000., 5000., 3000., 2000., 1000./
247	SAVE rlevstd
248	c$OMP THREADPRIVATE(rlevstd)
249	CHARACTER*4 clevSTD(nlevSTD)
250	DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
251	.'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
252	.'70 ','50 ','30 ','20 ','10 '/
253	SAVE clevSTD
254	c$OMP THREADPRIVATE(clevSTD)
255	c
256	CHARACTER*4 bb2
257	CHARACTER*2 bb3
258	c
259	real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD)
260	real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD)
261	real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)
262	real wlevSTD(klon,nlevSTD)
263
264	real twriteSTD(klon,nlevSTD,nfiles)
265	real qwriteSTD(klon,nlevSTD,nfiles)
266	real rhwriteSTD(klon,nlevSTD,nfiles)
267	real phiwriteSTD(klon,nlevSTD,nfiles)
268	real uwriteSTD(klon,nlevSTD,nfiles)
269	real vwriteSTD(klon,nlevSTD,nfiles)
270	real wwriteSTD(klon,nlevSTD,nfiles)
271	c
272	c nout : niveau de output des variables a une pression donnee
273	logical oknondef(klon,nlevSTD,nout)
274	c
275	c les produits uvSTD, vqSTD, .., T2STD sont calcules
276	c a partir des valeurs instantannees toutes les 6 h
277	c qui sont moyennees sur le mois
278	c
279	real uvSTD(klon,nlevSTD)
280	real vqSTD(klon,nlevSTD)
281	real vTSTD(klon,nlevSTD)
282	real wqSTD(klon,nlevSTD)
283	c
284	real vphiSTD(klon,nlevSTD)
285	real wTSTD(klon,nlevSTD)
286	real u2STD(klon,nlevSTD)
287	real v2STD(klon,nlevSTD)
288	real T2STD(klon,nlevSTD)
289	c
290	#include "radopt.h"
291	c
292	c
293	c prw: precipitable water
294	real prw(klon)
295
296	REAL convliq(klon,klev) ! eau liquide nuageuse convective
297	REAL convfra(klon,klev) ! fraction nuageuse convective
298
299	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
300	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
301	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
302	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
303
304	INTEGER linv, kp1
305	c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
306	c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
307	REAL flwp(klon), fiwp(klon)
308	REAL flwc(klon,klev), fiwc(klon,klev)
309	REAL flwp_c(klon), fiwp_c(klon)
310	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
311	REAL flwp_s(klon), fiwp_s(klon)
312	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
313
314	cIM ISCCP simulator v3.4
315	c dans clesphys.h top_height, overlap
316	cv3.4
317	INTEGER debug, debugcol
318	cym INTEGER npoints
319	cym PARAMETER(npoints=klon)
320	c
321	INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
322	INTEGER nregISCtot
323	PARAMETER(nregISCtot=1)
324	c
325	c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
326	c y compris pour 1 point
327	c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
328	c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
329	INTEGER imin_debut, nbpti
330	INTEGER jmin_debut, nbptj
331	cIM parametres ISCCP BEG
332	INTEGER nbapp_isccp
333	! INTEGER nbapp_isccp,isccppas
334	! PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
335	! !i.e. toutes les 3 heures
336	INTEGER n
337	INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
338	DATA ifreq_isccp/3/
339	SAVE ifreq_isccp
340	c$OMP THREADPRIVATE(ifreq_isccp)
341	CHARACTER*5 typinout(napisccp)
342	DATA typinout/'i3od'/
343	SAVE typinout
344	c$OMP THREADPRIVATE(typinout)
345	cIM verif boxptop BEG
346	CHARACTER*1 verticaxe(napisccp)
347	DATA verticaxe/'1'/
348	SAVE verticaxe
349	c$OMP THREADPRIVATE(verticaxe)
350	cIM verif boxptop END
351	INTEGER nvlev(napisccp)
352	c INTEGER nvlev
353	REAL t1, aa
354	REAL seed_re(klon,napisccp)
355	cym !!!! A voir plus tard
356	cym INTEGER iphy(iim,jjmp1)
357	cIM parametres ISCCP END
358	c
359	c ncol = nb. de sous-colonnes pour chaque maille du GCM
360	c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
361	c INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
362	INTEGER,SAVE :: ncol(napisccp)
363	c$OMP THREADPRIVATE(ncol)
364	INTEGER ncolmx, seed(klon,napisccp)
365	REAL nbsunlit(nregISCtot,klon,napisccp) !nbsunlit : moyenne de sunlit
366	c PARAMETER(ncolmx=1500)
367	PARAMETER(ncolmx=300)
368	c
369	cIM verif boxptop BEG
370	REAL vertlev(ncolmx,napisccp)
371	cIM verif boxptop END
372	c
373	REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
374	INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
375	c$OMP THREADPRIVATE(tautab,invtau)
376	REAL emsfc_lw
377	PARAMETER(emsfc_lw=0.99)
378	c REAL ran0 ! type for random number fuction
379	c
380	REAL cldtot(klon,klev)
381	c variables de haut en bas pour le simulateur ISCCP
382	REAL dtau_s(klon,klev) !tau nuages startiformes
383	REAL dtau_c(klon,klev) !tau nuages convectifs
384	REAL dem_s(klon,klev) !emissivite nuages startiformes
385	REAL dem_c(klon,klev) !emissivite nuages convectifs
386	c
387	c variables de haut en bas pour le simulateur ISCCP
388	REAL pfull(klon,klev)
389	REAL phalf(klon,klev+1)
390	REAL qv(klon,klev)
391	REAL cc(klon,klev)
392	REAL conv(klon,klev)
393	REAL dtau_sH2B(klon,klev)
394	REAL dtau_cH2B(klon,klev)
395	REAL at(klon,klev)
396	REAL dem_sH2B(klon,klev)
397	REAL dem_cH2B(klon,klev)
398	c
399	INTEGER kmax, lmax, lmax3
400	PARAMETER(kmax=8, lmax=8, lmax3=3)
401	INTEGER kmaxm1, lmaxm1
402	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
403	INTEGER iimx7, jjmx7, jjmp1x7
404	PARAMETER(iimx7=iimkmaxm1, jjmx7=jjmlmaxm1,
405	.jjmp1x7=jjmp1*lmaxm1)
406	c
407	c output from ISCCP simulator
408	REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
409	REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
410	REAL totalcldarea(klon,napisccp)
411	REAL meanptop(klon,napisccp)
412	REAL meantaucld(klon,napisccp)
413	REAL boxtau(klon,ncolmx,napisccp)
414	REAL boxptop(klon,ncolmx,napisccp)
415	REAL zx_tmp_fi3d_bx(klon,ncolmx)
416	REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
417	c
418	REAL cld_fi3d(klon,lmax3)
419	REAL cld_3d(iim,jjmp1,lmax3)
420	c
421	INTEGER iw, iwmax
422	REAL wmin, pas_w
423	c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
424	cIM 051005 PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
425	PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
426	REAL o500(klon)
427	c
428
429	c sorties ISCCP
430
431	integer nid_isccp
432	save nid_isccp
433	c$OMP THREADPRIVATE(nid_isccp)
434
435	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
436	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
437	SAVE zx_tau
438	DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
439	SAVE zx_pc
440	c$OMP THREADPRIVATE(zx_tau,zx_pc)
441	c cldtopres pression au sommet des nuages
442	REAL cldtopres(lmaxm1), cldtopres3(lmax3)
443	DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
444	DATA cldtopres3/440., 680., 1000./
445	SAVE cldtopres,cldtopres3
446	c$OMP THREADPRIVATE(cldtopres,cldtopres3)
447	cIM 051005 BEG
448	INTEGER komega, nhoriRD
449
450	c taulev: numero du niveau de tau dans les sorties ISCCP
451	CHARACTER *4 taulev(kmaxm1)
452	c DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
453	DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
454	CHARACTER *3 pclev(lmaxm1)
455	DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
456	SAVE taulev,pclev
457	c$OMP THREADPRIVATE(taulev,pclev)
458	c
459	c cnameisccp
460	CHARACTER *27 cnameisccp(lmaxm1,kmaxm1)
461	cIM bad 151205 DATA cnameisccp/'pc< 50hPa, tau< 0.3',
462	DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',
463	. 'pc= 180-310hPa, tau< 0.3',
464	. 'pc= 310-440hPa, tau< 0.3',
465	. 'pc= 440-560hPa, tau< 0.3',
466	. 'pc= 560-680hPa, tau< 0.3',
467	. 'pc= 680-800hPa, tau< 0.3',
468	. 'pc= 800-1000hPa, tau< 0.3',
469	. 'pc= 50-180hPa, tau= 0.3-1.3',
470	. 'pc= 180-310hPa, tau= 0.3-1.3',
471	. 'pc= 310-440hPa, tau= 0.3-1.3',
472	. 'pc= 440-560hPa, tau= 0.3-1.3',
473	. 'pc= 560-680hPa, tau= 0.3-1.3',
474	. 'pc= 680-800hPa, tau= 0.3-1.3',
475	. 'pc= 800-1000hPa, tau= 0.3-1.3',
476	. 'pc= 50-180hPa, tau= 1.3-3.6',
477	. 'pc= 180-310hPa, tau= 1.3-3.6',
478	. 'pc= 310-440hPa, tau= 1.3-3.6',
479	. 'pc= 440-560hPa, tau= 1.3-3.6',
480	. 'pc= 560-680hPa, tau= 1.3-3.6',
481	. 'pc= 680-800hPa, tau= 1.3-3.6',
482	. 'pc= 800-1000hPa, tau= 1.3-3.6',
483	. 'pc= 50-180hPa, tau= 3.6-9.4',
484	. 'pc= 180-310hPa, tau= 3.6-9.4',
485	. 'pc= 310-440hPa, tau= 3.6-9.4',
486	. 'pc= 440-560hPa, tau= 3.6-9.4',
487	. 'pc= 560-680hPa, tau= 3.6-9.4',
488	. 'pc= 680-800hPa, tau= 3.6-9.4',
489	. 'pc= 800-1000hPa, tau= 3.6-9.4',
490	. 'pc= 50-180hPa, tau= 9.4-23',
491	. 'pc= 180-310hPa, tau= 9.4-23',
492	. 'pc= 310-440hPa, tau= 9.4-23',
493	. 'pc= 440-560hPa, tau= 9.4-23',
494	. 'pc= 560-680hPa, tau= 9.4-23',
495	. 'pc= 680-800hPa, tau= 9.4-23',
496	. 'pc= 800-1000hPa, tau= 9.4-23',
497	. 'pc= 50-180hPa, tau= 23-60',
498	. 'pc= 180-310hPa, tau= 23-60',
499	. 'pc= 310-440hPa, tau= 23-60',
500	. 'pc= 440-560hPa, tau= 23-60',
501	. 'pc= 560-680hPa, tau= 23-60',
502	. 'pc= 680-800hPa, tau= 23-60',
503	. 'pc= 800-1000hPa, tau= 23-60',
504	. 'pc= 50-180hPa, tau> 60.',
505	. 'pc= 180-310hPa, tau> 60.',
506	. 'pc= 310-440hPa, tau> 60.',
507	. 'pc= 440-560hPa, tau> 60.',
508	. 'pc= 560-680hPa, tau> 60.',
509	. 'pc= 680-800hPa, tau> 60.',
510	. 'pc= 800-1000hPa, tau> 60.'/
511	SAVE cnameisccp
512	c$OMP THREADPRIVATE(cnameisccp)
513	c
514	c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
515	c INTEGER nhorix7
516	cIM: region='3d' <==> sorties en global
517	CHARACTER*3 region
518	PARAMETER(region='3d')
519	c
520	cIM ISCCP simulator v3.4
521	c
522	logical ok_hf
523	c
524	integer nid_hf, nid_hf3d
525	save ok_hf, nid_hf, nid_hf3d
526	c$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
527	c QUESTION : noms de variables ?
528
529	INTEGER longcles
530	PARAMETER ( longcles = 20 )
531	REAL clesphy0( longcles )
532	c
533	c Variables propres a la physique
534	INTEGER itap
535	SAVE itap ! compteur pour la physique
536	c$OMP THREADPRIVATE(itap)
537	c
538	real slp(klon) ! sea level pressure
539	c
540	REAL fevap(klon,nbsrf)
541	REAL fluxlat(klon,nbsrf)
542	c
543	REAL qsol(klon)
544	REAL,save :: solarlong0
545	c$OMP THREADPRIVATE(solarlong0)
546
547	c
548	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
549	c
550	cIM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
551	REAL zulow(klon),zvlow(klon)
552	c
553	INTEGER igwd,idx(klon),itest(klon)
554	c
555	REAL agesno(klon,nbsrf)
556	c
557	c REAL,allocatable,save :: run_off_lic_0(:)
558	cc$OMP THREADPRIVATE(run_off_lic_0)
559	cym SAVE run_off_lic_0
560	cKE43
561	c Variables liees a la convection de K. Emanuel (sb):
562	c
563	REAL bas, top ! cloud base and top levels
564	SAVE bas
565	SAVE top
566	c$OMP THREADPRIVATE(bas, top)
567
568	REAL wdn(klon), tdn(klon), qdn(klon)
569	c
570	c=================================================================================================
571	cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
572	c Variables liées à la poche froide (jyg)
573
574	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
575	REAL Vprecip(klon,klev) ! precipitation vertical profile
576	c
577	REAL wape_prescr, fip_prescr
578	INTEGER it_wape_prescr
579	SAVE wape_prescr, fip_prescr, it_wape_prescr
580	c$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
581	c
582	c variables supplementaires de concvl
583	REAL Tconv(klon,klev)
584	REAL ment(klon,klev,klev),sij(klon,klev,klev)
585	REAL dd_t(klon,klev),dd_q(klon,klev)
586
587	real, save :: alp_bl_prescr=0.
588	real, save :: ale_bl_prescr=0.
589
590	real, save :: ale_max=100.
591	real, save :: alp_max=2.
592
593	c$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
594	c$OMP THREADPRIVATE(ale_max,alp_max)
595
596	real ale_wake(klon)
597	real alp_wake(klon)
598	cRC
599	c Variables liées à la poche froide (jyg et rr)
600	c Version diagnostique pour l'instant : pas de rétroaction sur la convection
601
602	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
603
604	REAL wake_dth(klon,klev) ! wake : temp pot difference
605
606	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
607	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
608	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
609	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
610	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
611	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
612	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
613	REAL wake_omg(klon,klev) ! Wake : velocity difference (wake - unpertubed)
614	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
615	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
616	REAL wake_spread(klon,klev) ! spreading term in wake_delt
617	c
618	cpourquoi y'a pas de save??
619	REAL wake_h(klon) ! Wake : hauteur de la poche froide
620	c
621	INTEGER wake_k(klon) ! Wake sommet
622	c
623	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
624	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
625	c
626	REAL wake_pe(klon) ! Wake potential energy - WAPE
627
628	REAL wake_gfl(klon) ! Gust Front Length
629	REAL wake_dens(klon)
630	c
631	c
632	REAL dt_dwn(klon,klev)
633	REAL dq_dwn(klon,klev)
634	REAL wdt_PBL(klon,klev)
635	REAL udt_PBL(klon,klev)
636	REAL wdq_PBL(klon,klev)
637	REAL udq_PBL(klon,klev)
638	REAL M_dwn(klon,klev)
639	REAL M_up(klon,klev)
640	REAL dt_a(klon,klev)
641	REAL dq_a(klon,klev)
642	c
643	cRR:fin declarations poches froides
644	c=======================================================================================================
645
646	REAL zw2(klon,klev+1)
647	REAL fraca(klon,klev+1)
648
649	c Variables locales pour la couche limite (al1):
650	c
651	cAl1 REAL pblh(klon) ! Hauteur de couche limite
652	cAl1 SAVE pblh
653	c34EK
654	c
655	c Variables locales:
656	c
657	REAL cdragh(klon) ! drag coefficient pour T and Q
658	REAL cdragm(klon) ! drag coefficient pour vent
659	cAA
660	cAA Pour phytrac
661	cAA
662	REAL coefh(klon,klev) ! coef d'echange pour phytrac, valable pour 2<=k<=klev
663	REAL u1(klon) ! vents dans la premiere couche U
664	REAL v1(klon) ! vents dans la premiere couche V
665
666	REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
667
668	c@$$ LOGICAL offline ! Controle du stockage ds "physique"
669	c@$$ PARAMETER (offline=.false.)
670	c@$$ INTEGER physid
671	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
672	REAL frac_nucl(klon,klev) ! idem (nucleation)
673	INTEGER :: iii
674	REAL :: calday
675
676	cIM cf FH pour Tiedtke 080604
677	REAL rain_tiedtke(klon),snow_tiedtke(klon)
678	c
679	cIM 050204 END
680	REAL evap(klon), devap(klon) ! evaporation et sa derivee
681	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
682
683	REAL bils(klon) ! bilan de chaleur au sol
684	REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
685	C ! type de sous-surface et pondere par la fraction
686	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
687	C ! type de sous-surface et pondere par la fraction
688	REAL slab_wfbils(klon) ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
689
690	REAL fder(klon)
691	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
692	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
693	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
694	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
695	c
696	REAL frugs(klon,nbsrf)
697	REAL zxrugs(klon) ! longueur de rugosite
698	c
699	c Conditions aux limites
700	c
701	!
702	REAL :: day_since_equinox
703	! Date de l'equinoxe de printemps
704	INTEGER, parameter :: mth_eq=3, day_eq=21
705	REAL :: jD_eq
706
707	LOGICAL, parameter :: new_orbit = .true.
708
709	c
710	INTEGER lmt_pas
711	SAVE lmt_pas ! frequence de mise a jour
712	c$OMP THREADPRIVATE(lmt_pas)
713	real zmasse(klon, llm)
714	C (column-density of mass of air in a cell, in kg m-2)
715	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
716
717	cIM sorties
718	REAL un_jour
719	PARAMETER(un_jour=86400.)
720	c======================================================================
721	c
722	c Declaration des procedures appelees
723	c
724	EXTERNAL angle ! calculer angle zenithal du soleil
725	EXTERNAL alboc ! calculer l'albedo sur ocean
726	EXTERNAL ajsec ! ajustement sec
727	EXTERNAL conlmd ! convection (schema LMD)
728	cKE43
729	EXTERNAL conema3 ! convect4.3
730	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
731	cAA
732	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
733	c ! stockage des coefficients necessaires au
734	c ! lessivage OFF-LINE et ON-LINE
735	EXTERNAL hgardfou ! verifier les temperatures
736	EXTERNAL nuage ! calculer les proprietes radiatives
737	CC EXTERNAL o3cm ! initialiser l'ozone
738	EXTERNAL orbite ! calculer l'orbite terrestre
739	EXTERNAL phyetat0 ! lire l'etat initial de la physique
740	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
741	EXTERNAL suphel ! initialiser certaines constantes
742	EXTERNAL transp ! transport total de l'eau et de l'energie
743	EXTERNAL ecribina ! ecrire le fichier binaire global
744	EXTERNAL ecribins ! ecrire le fichier binaire global
745	EXTERNAL ecrirega ! ecrire le fichier binaire regional
746	EXTERNAL ecriregs ! ecrire le fichier binaire regional
747	cIM
748	EXTERNAL haut2bas !variables de haut en bas
749	INTEGER lnblnk1
750	EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type
751	!caracter
752	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
753	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
754	c EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
755	c EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
756	c !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
757	c
758	c Variables locales
759	c
760	REAL rhcl(klon,klev) ! humiditi relative ciel clair
761	REAL dialiq(klon,klev) ! eau liquide nuageuse
762	REAL diafra(klon,klev) ! fraction nuageuse
763	REAL cldliq(klon,klev) ! eau liquide nuageuse
764	REAL cldfra(klon,klev) ! fraction nuageuse
765	REAL cldtau(klon,klev) ! epaisseur optique
766	REAL cldemi(klon,klev) ! emissivite infrarouge
767	c
768	CXXX PB
769	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
770	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
771	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
772	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
773	c
774	REAL zxfluxt(klon, klev)
775	REAL zxfluxq(klon, klev)
776	REAL zxfluxu(klon, klev)
777	REAL zxfluxv(klon, klev)
778	CXXX
779	c
780	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
781	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
782	c Le rayonnement n'est pas calcule tous les pas, il faut donc
783	c sauvegarder les sorties du rayonnement
784	cym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
785	cym SAVE sollwdownclr, toplwdown, toplwdownclr
786	cym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
787	c
788	INTEGER itaprad
789	SAVE itaprad
790	c$OMP THREADPRIVATE(itaprad)
791	c
792	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
793	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
794	c
795	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
796	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
797	c
798	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
799	REAL zxsnow_dummy(klon)
800	c
801	REAL dist, rmu0(klon), fract(klon)
802	REAL zdtime, zlongi
803	c
804	CHARACTER*2 str2
805	CHARACTER*2 iqn
806	c
807	REAL qcheck
808	REAL z_avant(klon), z_apres(klon), z_factor(klon)
809	LOGICAL zx_ajustq
810	c
811	REAL za, zb
812	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
813	real zqsat(klon,klev)
814	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
815	REAL t_coup
816	PARAMETER (t_coup=234.0)
817	c
818	REAL zphi(klon,klev)
819	cym A voir plus tard !!
820	cym REAL zx_relief(iim,jjmp1)
821	cym REAL zx_aire(iim,jjmp1)
822	c
823	c Grandeurs de sorties
824	REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
825	REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
826	REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
827	REAL s_trmb3(klon)
828	cKE43
829	c Variables locales pour la convection de K. Emanuel (sb):
830	c
831	REAL upwd(klon,klev) ! saturated updraft mass flux
832	REAL dnwd(klon,klev) ! saturated downdraft mass flux
833	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
834	REAL tvp(klon,klev) ! virtual temp of lifted parcel
835	CHARACTER*40 capemaxcels !max(CAPE)
836
837	REAL rflag(klon) ! flag fonctionnement de convect
838	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
839	c -- convect43:
840	INTEGER ntra ! nb traceurs pour convect4.3
841	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
842	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
843	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
844	REAL dplcldt(klon), dplcldr(klon)
845	c? . condm_con(klon,klev),conda_con(klon,klev),
846	c? . mr_con(klon,klev),ep_con(klon,klev)
847	c? . ,sadiab(klon,klev),wadiab(klon,klev)
848	c --
849	c34EK
850	c
851	c Variables du changement
852	c
853	c con: convection
854	c lsc: condensation a grande echelle (Large-Scale-Condensation)
855	c ajs: ajustement sec
856	c eva: evaporation de l'eau liquide nuageuse
857	c vdf: couche limite (Vertical DiFfusion)
858	REAL rneb(klon,klev)
859
860	! tendance nulles
861	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
862
863	c
864	*********************************************************
865	* declarations
866
867	*********************************************************
868	cIM 081204 END
869	c
870	REAL pmfu(klon,klev), pmfd(klon,klev)
871	REAL pen_u(klon,klev), pen_d(klon,klev)
872	REAL pde_u(klon,klev), pde_d(klon,klev)
873	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
874	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
875	REAL prfl(klon,klev+1), psfl(klon,klev+1)
876	c
877	REAL rain_lsc(klon)
878	REAL snow_lsc(klon)
879	c
880	REAL ratqss(klon,klev),ratqsc(klon,klev)
881	real ratqsbas,ratqshaut,tau_ratqs
882	save ratqsbas,ratqshaut,tau_ratqs
883	c$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
884	real zpt_conv(klon,klev)
885
886	c Parametres lies au nouveau schema de nuages (SB, PDF)
887	real fact_cldcon
888	real facttemps
889	logical ok_newmicro
890	save ok_newmicro
891	real ref_liq(klon,klev), ref_ice(klon,klev)
892	c$OMP THREADPRIVATE(ok_newmicro)
893	save fact_cldcon,facttemps
894	c$OMP THREADPRIVATE(fact_cldcon,facttemps)
895	real facteur
896
897	integer iflag_cldcon
898	save iflag_cldcon
899	c$OMP THREADPRIVATE(iflag_cldcon)
900	logical ptconv(klon,klev)
901	cIM cf. AM 081204 BEG
902	logical ptconvth(klon,klev)
903	cIM cf. AM 081204 END
904	c
905	c Variables liees a l'ecriture de la bande histoire physique
906	c
907	c======================================================================
908	c
909	cIM cf. AM 081204 BEG
910	c declarations pour sortir sur une sous-region
911	integer imin_ins,imax_ins,jmin_ins,jmax_ins
912	save imin_ins,imax_ins,jmin_ins,jmax_ins
913	c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
914	c real lonmin_ins,lonmax_ins,latmin_ins
915	c s ,latmax_ins
916	c data lonmin_ins,lonmax_ins,latmin_ins
917	c s ,latmax_ins/
918	c valeurs initiales s -5.,20.,41.,55./
919	c s 100.,130.,-20.,20./
920	c s -180.,180.,-90.,90./
921	c======================================================================
922	cIM cf. AM 081204 END
923
924	c
925	integer itau_w ! pas de temps ecriture = itap + itau_phy
926	c
927	c
928	c Variables locales pour effectuer les appels en serie
929	c
930	REAL zx_rh(klon,klev)
931	cIM RH a 2m (la surface)
932	REAL rh2m(klon), qsat2m(klon)
933	REAL tpot(klon), tpote(klon)
934	REAL Lheat
935
936	INTEGER length
937	PARAMETER ( length = 100 )
938	REAL tabcntr0( length )
939	c
940	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
941	cIM
942	INTEGER ndex2d1(iwmax)
943	c
944	cIM AMIP2 BEG
945	REAL moyglo, mountor
946	cIM 141004 BEG
947	REAL zustrdr(klon), zvstrdr(klon)
948	REAL zustrli(klon), zvstrli(klon)
949	REAL zustrph(klon), zvstrph(klon)
950	REAL zustrhi(klon), zvstrhi(klon)
951	REAL aam, torsfc
952	cIM 141004 END
953	cIM 190504 BEG
954	INTEGER ij, imp1jmp1
955	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
956	cym A voir plus tard
957	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
958	REAL padyn(iim+1,jjmp1,klev+1)
959	REAL dudyn(iim+1,jjmp1,klev)
960	REAL rlatdyn(iim+1,jjmp1)
961	cIM 190504 END
962	LOGICAL ok_msk
963	REAL msk(klon)
964	cIM
965	REAL airetot, pi
966	cym A voir plus tard
967	cym REAL zm_wo(jjmp1, klev)
968	cIM AMIP2 END
969	c
970	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
971	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
972	c#ifdef histmthNMC
973	cym A voir plus tard !!!!
974	cym REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
975	REAL zx_tmp_fiNC(klon,nlevSTD)
976	c#endif
977	REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
978	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
979	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
980	c
981	INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
982	INTEGER nid_ctesGCM
983	SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
984	SAVE nid_ctesGCM
985	c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins, nid_nmc)
986	c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
987	c
988	cIM 280405 BEG
989	INTEGER nid_bilKPins, nid_bilKPave
990	SAVE nid_bilKPins, nid_bilKPave
991	c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
992	c
993	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
994	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
995	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
996	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
997	c
998	cIM 280405 END
999	c
1000	INTEGER nhori, nvert, nvert1, nvert3
1001	REAL zsto, zsto1, zsto2
1002	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1003	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1004	REAL zout_isccp(napisccp)
1005	SAVE zcals, zcalh, zoutj, zout_isccp
1006	c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1007
1008	real zjulian
1009	save zjulian
1010	c$OMP THREADPRIVATE(zjulian)
1011
1012	character*20 modname
1013	character*80 abort_message
1014	logical ok_sync
1015	real date0
1016	integer idayref
1017
1018	C essai writephys
1019	integer fid_day, fid_mth, fid_ins
1020	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1021	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1022	parameter (prof2d_on = 1, prof3d_on = 2,
1023	. prof2d_av = 3, prof3d_av = 4)
1024	character*30 nom_fichier
1025	character*10 varname
1026	character*40 vartitle
1027	character*20 varunits
1028	C Variables liees au bilan d'energie et d'enthalpi
1029	REAL ztsol(klon)
1030	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1031	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1032	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1033	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1034	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1035	c$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1036	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1037	REAL d_h_vcol_phy
1038	REAL fs_bound, fq_bound
1039	SAVE d_h_vcol_phy
1040	c$OMP THREADPRIVATE(d_h_vcol_phy)
1041	REAL zero_v(klon)
1042	CHARACTER*15 ztit
1043	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
1044	SAVE ip_ebil
1045	DATA ip_ebil/0/
1046	c$OMP THREADPRIVATE(ip_ebil)
1047	INTEGER if_ebil ! level for energy conserv. dignostics
1048	SAVE if_ebil
1049	c$OMP THREADPRIVATE(if_ebil)
1050	c+jld ec_conser
1051	REAL ZRCPD
1052	c-jld ec_conser
1053	REAL t2m(klon,nbsrf) ! temperature a 2m
1054	REAL q2m(klon,nbsrf) ! humidite a 2m
1055
1056	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1057	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1058	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
1059	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1060	CHARACTER*40 tinst, tave, typeval
1061	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1062
1063	REAL re(klon, klev) ! Cloud droplet effective radius
1064	REAL fl(klon, klev) ! denominator of re
1065
1066	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1067
1068	! Aerosol optical properties
1069	CHARACTER*4, DIMENSION(naero_grp) :: rfname
1070	REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index
1071	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass concentration for all soluble aerosols[ug/m3]
1072	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi ! - " - (pre-industrial value)
1073	INTEGER :: naero ! aerosol species
1074
1075	! Parameters
1076	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1077	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1078	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1079	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
1080	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1081	! false : lecture des aerosol dans un fichier
1082	c$OMP THREADPRIVATE(aerosol_couple)
1083	INTEGER, SAVE :: flag_aerosol
1084	c$OMP THREADPRIVATE(flag_aerosol)
1085	LOGICAL, SAVE :: new_aod
1086	c$OMP THREADPRIVATE(new_aod)
1087
1088	c
1089	c Declaration des constantes et des fonctions thermodynamiques
1090	c
1091	LOGICAL,SAVE :: first=.true.
1092	c$OMP THREADPRIVATE(first)
1093
1094	integer iunit
1095
1096	integer, save:: read_climoz ! read ozone climatology
1097	C Allowed values are 0, 1 and 2
1098	C 0: do not read an ozone climatology
1099	C 1: read a single ozone climatology that will be used day and night
1100	C 2: read two ozone climatologies, the average day and night
1101	C climatology and the daylight climatology
1102
1103	integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
1104
1105	real, pointer, save:: press_climoz(:)
1106	! edges of pressure intervals for ozone climatologies, in Pa, in strictly
1107	! ascending order
1108
1109	integer, save:: co3i = 0
1110	! time index in NetCDF file of current ozone fields
1111	c$OMP THREADPRIVATE(co3i)
1112
1113	integer ro3i
1114	! required time index in NetCDF file for the ozone fields, between 1
1115	! and 360
1116
1117	#include "YOMCST.h"
1118	#include "YOETHF.h"
1119	#include "FCTTRE.h"
1120	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1121	#include "conema3.h"
1122	#include "fisrtilp.h"
1123	#include "nuage.h"
1124	#include "compbl.h"
1125	cIM 100106 END : pouvoir sortir les ctes de la physique
1126	c
1127	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1128	c Declarations pour Simulateur COSP
1129	c============================================================
1130	real :: mr_ozone(klon,klev)
1131	c======================================================================
1132	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1133	! Utilise notamment en 1D mais peut etre active egalement en 3D
1134	! en imposant la valeur de igout.
1135	c======================================================================
1136
1137	if (prt_level.ge.1) then
1138	igout=klon/2+1/klon
1139	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1140	write(lunout,*)
1141	s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1142	write(lunout,*)
1143	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1144
1145	write(lunout,*) 'paprs, play, phi, u, v, t'
1146	do k=1,klev
1147	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1148	s u(igout,k),v(igout,k),t(igout,k)
1149	enddo
1150	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1151	do k=1,klev
1152	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1153	enddo
1154	endif
1155
1156	c======================================================================
1157
1158	cym => necessaire pour iflag_con != 2
1159	pmfd(:,:) = 0.
1160	pen_u(:,:) = 0.
1161	pen_d(:,:) = 0.
1162	pde_d(:,:) = 0.
1163	pde_u(:,:) = 0.
1164	aam=0.
1165
1166	torsfc=0.
1167	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1168
1169	if (first) then
1170
1171	cCR:nvelles variables convection/poches froides
1172
1173	print*, '================================================='
1174	print*, 'Allocation des variables locales et sauvegardees'
1175	call phys_local_var_init
1176	c appel a la lecture du run.def physique
1177	call conf_phys(ok_journe, ok_mensuel,
1178	. ok_instan, ok_hf,
1179	. ok_LES,
1180	. solarlong0,seuil_inversion,
1181	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1182	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
1183	. ok_ade, ok_aie, aerosol_couple,
1184	. flag_aerosol, new_aod,
1185	. bl95_b0, bl95_b1,
1186	. iflag_thermals,nsplit_thermals,tau_thermals,
1187	. iflag_thermals_ed,iflag_thermals_optflux,
1188	c nv flags pour la convection et les poches froides
1189	. iflag_coupl,iflag_clos,iflag_wake, read_climoz)
1190	call phys_state_var_init(read_climoz)
1191	print*, '================================================='
1192
1193	cIM beg
1194	dnwd0=0.0
1195	ftd=0.0
1196	fqd=0.0
1197	cin=0.
1198	cym Attention pbase pas initialise dans concvl !!!!
1199	pbase=0
1200	paire_ter(:)=0.
1201	cIM 180608
1202	c pmflxr=0.
1203	c pmflxs=0.
1204	first=.false.
1205
1206	endif ! first
1207
1208	modname = 'physiq'
1209	cIM
1210	IF (ip_ebil_phy.ge.1) THEN
1211	DO i=1,klon
1212	zero_v(i)=0.
1213	END DO
1214	END IF
1215	ok_sync=.TRUE.
1216
1217	IF (debut) THEN
1218	CALL suphel ! initialiser constantes et parametres phys.
1219	ENDIF
1220
1221	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1222
1223
1224	c======================================================================
1225	! Gestion calendrier : mise a jour du module phys_cal_mod
1226	!
1227	CALL phys_cal_update(jD_cur,jH_cur)
1228
1229	c
1230	c Si c'est le debut, il faut initialiser plusieurs choses
1231	c ********
1232	c
1233	IF (debut) THEN
1234	!rv
1235	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1236	cde la convection a partir des caracteristiques du thermique
1237	wght_th(:,:)=1.
1238	lalim_conv(:)=1
1239	cRC
1240	u10m(:,:)=0.
1241	v10m(:,:)=0.
1242	rain_con(:)=0.
1243	snow_con(:)=0.
1244	topswai(:)=0.
1245	topswad(:)=0.
1246	solswai(:)=0.
1247	solswad(:)=0.
1248
1249	lambda_th(:,:)=0.
1250	wmax_th(:)=0.
1251	tau_overturning_th(:)=0.
1252
1253	IF (config_inca /= 'none') THEN
1254	! jg : initialisation jusqu'au ces variables sont dans restart
1255	ccm(:,:,:) = 0.
1256	tau_aero(:,:,:,:) = 0.
1257	piz_aero(:,:,:,:) = 0.
1258	cg_aero(:,:,:,:) = 0.
1259	END IF
1260
1261	rnebcon0(:,:) = 0.0
1262	clwcon0(:,:) = 0.0
1263	rnebcon(:,:) = 0.0
1264	clwcon(:,:) = 0.0
1265
1266	cIM
1267	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1268	c
1269	print*,'iflag_coupl,iflag_clos,iflag_wake',
1270	. iflag_coupl,iflag_clos,iflag_wake
1271	print*,'CYCLE_DIURNE', cycle_diurne
1272	c
1273	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1274	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1275	CALL abort_gcm (modname,abort_message,1)
1276	ENDIF
1277	c
1278	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1279	abort_message = 'ISCCP-like outputs may be available for KE
1280	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1281	CALL abort_gcm (modname,abort_message,1)
1282	ENDIF
1283	c
1284	c Initialiser les compteurs:
1285	c
1286	itap = 0
1287	itaprad = 0
1288
1289	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1290	!! Un petit travail à faire ici.
1291	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1292
1293	if (iflag_pbl>1) then
1294	PRINT*, "Using method MELLOR&YAMADA"
1295	endif
1296
1297	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1298	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1299	! dyn3d
1300	! Attention : la version precedente n'etait pas tres propre.
1301	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1302	! pour obtenir le meme resultat.
1303	dtime=pdtphys
1304	radpas = NINT( 86400./dtime/nbapp_rad)
1305	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1306
1307	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1308	cIM begin
1309	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1310	$,ratqs(1,1)
1311	cIM end
1312
1313
1314
1315	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1316	c
1317	C on remet le calendrier a zero
1318	c
1319	IF (raz_date .eq. 1) THEN
1320	itau_phy = 0
1321	ENDIF
1322
1323	cIM cf. AM 081204 BEG
1324	PRINT*,'cycle_diurne3 =',cycle_diurne
1325	cIM cf. AM 081204 END
1326	c
1327	CALL printflag( tabcntr0,radpas,ok_journe,
1328	, ok_instan, ok_region )
1329	c
1330	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1331	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1332	. pdtphys
1333	abort_message='Pas physique n est pas correct '
1334	! call abort_gcm(modname,abort_message,1)
1335	dtime=pdtphys
1336	ENDIF
1337	IF (nlon .NE. klon) THEN
1338	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1339	. klon
1340	abort_message='nlon et klon ne sont pas coherents'
1341	call abort_gcm(modname,abort_message,1)
1342	ENDIF
1343	IF (nlev .NE. klev) THEN
1344	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1345	. klev
1346	abort_message='nlev et klev ne sont pas coherents'
1347	call abort_gcm(modname,abort_message,1)
1348	ENDIF
1349	c
1350	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1351	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1352	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1353	abort_message='Nbre d appels au rayonnement insuffisant'
1354	call abort_gcm(modname,abort_message,1)
1355	ENDIF
1356	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1357	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1358	. ok_cvl
1359	c
1360	cKE43
1361	c Initialisation pour la convection de K.E. (sb):
1362	IF (iflag_con.GE.3) THEN
1363
1364	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1365	WRITE(lunout,*)
1366	. "On va utiliser le melange convectif des traceurs qui"
1367	WRITE(lunout,*)"est calcule dans convect4.3"
1368	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1369
1370	DO i = 1, klon
1371	ema_cbmf(i) = 0.
1372	ema_pcb(i) = 0.
1373	ema_pct(i) = 0.
1374	ema_workcbmf(i) = 0.
1375	ENDDO
1376	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1377	DO i = 1, klon
1378	ibas_con(i) = 1
1379	itop_con(i) = 1
1380	ENDDO
1381	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1382	c===============================================================================
1383	cCR:04.12.07: initialisations poches froides
1384	c Controle de ALE et ALP pour la fermeture convective (jyg)
1385	if (iflag_wake.eq.1) then
1386	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1387	s ,alp_bl_prescr, ale_bl_prescr)
1388	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1389	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1390	endif
1391
1392	do i = 1,klon
1393	Ale_bl(i)=0.
1394	Alp_bl(i)=0.
1395	enddo
1396
1397	c================================================================================
1398
1399	ENDIF
1400
1401	DO i=1,klon
1402	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1403	ENDDO
1404
1405	c34EK
1406	IF (ok_orodr) THEN
1407
1408	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1409	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1410	! justement quand ok_orodr = false.
1411	! ce rugoro est utilise par la couche limite et fait double emploi
1412	! avec les paramétrisations spécifiques de Francois Lott.
1413	! DO i=1,klon
1414	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1415	! ENDDO
1416	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1417	IF (ok_strato) THEN
1418	CALL SUGWD_strato(klon,klev,paprs,pplay)
1419	ELSE
1420	CALL SUGWD(klon,klev,paprs,pplay)
1421	ENDIF
1422
1423	DO i=1,klon
1424	zuthe(i)=0.
1425	zvthe(i)=0.
1426	if(zstd(i).gt.10.)then
1427	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1428	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1429	endif
1430	ENDDO
1431	ENDIF
1432	c
1433	c
1434	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1435	WRITE(lunout,*)'La frequence de lecture surface est de ',
1436	. lmt_pas
1437	c
1438	cIM 030306 END
1439
1440	capemaxcels = 't_max(X)'
1441	t2mincels = 't_min(X)'
1442	t2maxcels = 't_max(X)'
1443	tinst = 'inst(X)'
1444	tave = 'ave(X)'
1445	cIM cf. AM 081204 BEG
1446	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1447	cIM cf. AM 081204 END
1448	c
1449	c=============================================================
1450	c Initialisation des sorties
1451	c=============================================================
1452
1453	#ifdef CPP_IOIPSL
1454
1455	c$OMP MASTER
1456	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1457	& ctetaSTD,dtime,ok_veget,
1458	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1459	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, read_climoz)
1460	c$OMP END MASTER
1461	c$OMP BARRIER
1462
1463	#ifdef histISCCP
1464	#include "ini_histISCCP.h"
1465	#endif
1466
1467	#ifdef histmthNMC
1468	#include "ini_histmthNMC.h"
1469	#endif
1470
1471	#include "ini_histday_seri.h"
1472
1473	#include "ini_paramLMDZ_phy.h"
1474
1475	#endif
1476
1477	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1478	ecrit_hf2mth = ecrit_mth/ecrit_hf
1479
1480	ecrit_hf = ecrit_hf * un_jour
1481	!IM
1482	IF(ecrit_day.LE.1.) THEN
1483	ecrit_day = ecrit_day * un_jour !en secondes
1484	ENDIF
1485	!IM
1486	ecrit_mth = ecrit_mth * un_jour
1487	ecrit_ins = ecrit_ins * un_jour
1488	ecrit_reg = ecrit_reg * un_jour
1489	ecrit_tra = ecrit_tra * un_jour
1490	ecrit_ISCCP = ecrit_ISCCP * un_jour
1491	ecrit_LES = ecrit_LES * un_jour
1492	c
1493	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1494	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1495	. ecrit_hf2mth
1496	cIM 030306 END
1497
1498
1499	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1500	date0 = jD_ref
1501	WRITE(,) 'physiq date0 : ',date0
1502	c
1503	c
1504	c
1505	c Prescrire l'ozone dans l'atmosphere
1506	c
1507	c
1508	cc DO i = 1, klon
1509	cc DO k = 1, klev
1510	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1511	cc ENDDO
1512	cc ENDDO
1513	c
1514	IF (config_inca /= 'none') THEN
1515	#ifdef INCA
1516	CALL VTe(VTphysiq)
1517	CALL VTb(VTinca)
1518	! iii = MOD(NINT(xjour),360)
1519	! calday = FLOAT(iii) + jH_cur
1520	calday = FLOAT(days_elapsed) + jH_cur
1521	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1522
1523	CALL chemini(
1524	$ rg,
1525	$ ra,
1526	$ airephy,
1527	$ rlat,
1528	$ rlon,
1529	$ presnivs,
1530	$ calday,
1531	$ klon,
1532	$ nqtot,
1533	$ pdtphys,
1534	$ annee_ref,
1535	$ day_ini)
1536
1537	CALL VTe(VTinca)
1538	CALL VTb(VTphysiq)
1539	#endif
1540	END IF
1541	c
1542	c
1543	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1544	! Nouvelle initialisation pour le rayonnement RRTM
1545	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1546
1547	call iniradia(klon,klev,paprs(1,1:klev+1))
1548
1549	C$omp single
1550	if (read_climoz >= 1) then
1551	call open_climoz(ncid_climoz, press_climoz)
1552	END IF
1553	C$omp end single
1554	ENDIF
1555	!
1556	! ************** Fin de IF ( debut ) *************
1557	!
1558	!
1559	! Incrementer le compteur de la physique
1560	!
1561	itap = itap + 1
1562	!
1563	! Update fraction of the sub-surfaces (pctsrf) and
1564	! initialize, where a new fraction has appeared, all variables depending
1565	! on the surface fraction.
1566	!
1567	CALL change_srf_frac(itap, dtime, days_elapsed+1,
1568	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1569
1570	! Tendances bidons pour les processus qui n'affectent pas certaines
1571	! variables.
1572	du0(:,:)=0.
1573	dv0(:,:)=0.
1574	dq0(:,:)=0.
1575	dql0(:,:)=0.
1576	c
1577	c Mettre a zero des variables de sortie (pour securite)
1578	c
1579	DO i = 1, klon
1580	d_ps(i) = 0.0
1581	ENDDO
1582	DO k = 1, klev
1583	DO i = 1, klon
1584	d_t(i,k) = 0.0
1585	d_u(i,k) = 0.0
1586	d_v(i,k) = 0.0
1587	ENDDO
1588	ENDDO
1589	DO iq = 1, nqtot
1590	DO k = 1, klev
1591	DO i = 1, klon
1592	d_qx(i,k,iq) = 0.0
1593	ENDDO
1594	ENDDO
1595	ENDDO
1596	da(:,:)=0.
1597	mp(:,:)=0.
1598	phi(:,:,:)=0.
1599	c
1600	c Ne pas affecter les valeurs entrees de u, v, h, et q
1601	c
1602	DO k = 1, klev
1603	DO i = 1, klon
1604	t_seri(i,k) = t(i,k)
1605	u_seri(i,k) = u(i,k)
1606	v_seri(i,k) = v(i,k)
1607	q_seri(i,k) = qx(i,k,ivap)
1608	ql_seri(i,k) = qx(i,k,iliq)
1609	qs_seri(i,k) = 0.
1610	ENDDO
1611	ENDDO
1612	IF (nqtot.GE.3) THEN
1613	DO iq = 3, nqtot
1614	DO k = 1, klev
1615	DO i = 1, klon
1616	tr_seri(i,k,iq-2) = qx(i,k,iq)
1617	ENDDO
1618	ENDDO
1619	ENDDO
1620	ELSE
1621	DO k = 1, klev
1622	DO i = 1, klon
1623	tr_seri(i,k,1) = 0.0
1624	ENDDO
1625	ENDDO
1626	ENDIF
1627	C
1628	DO i = 1, klon
1629	ztsol(i) = 0.
1630	ENDDO
1631	DO nsrf = 1, nbsrf
1632	DO i = 1, klon
1633	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1634	ENDDO
1635	ENDDO
1636	cIM
1637	IF (ip_ebil_phy.ge.1) THEN
1638	ztit='after dynamic'
1639	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1640	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1641	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1642	C Comme les tendances de la physique sont ajoute dans la dynamique,
1643	C on devrait avoir que la variation d'entalpie par la dynamique
1644	C est egale a la variation de la physique au pas de temps precedent.
1645	C Donc la somme de ces 2 variations devrait etre nulle.
1646	call diagphy(airephy,ztit,ip_ebil_phy
1647	e , zero_v, zero_v, zero_v, zero_v, zero_v
1648	e , zero_v, zero_v, zero_v, ztsol
1649	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1650	s , fs_bound, fq_bound )
1651	END IF
1652
1653	c Diagnostiquer la tendance dynamique
1654	c
1655	IF (ancien_ok) THEN
1656	DO k = 1, klev
1657	DO i = 1, klon
1658	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1659	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1660	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1661	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1662	ENDDO
1663	ENDDO
1664	ELSE
1665	DO k = 1, klev
1666	DO i = 1, klon
1667	d_u_dyn(i,k) = 0.0
1668	d_v_dyn(i,k) = 0.0
1669	d_t_dyn(i,k) = 0.0
1670	d_q_dyn(i,k) = 0.0
1671	ENDDO
1672	ENDDO
1673	ancien_ok = .TRUE.
1674	ENDIF
1675	c
1676	c Ajouter le geopotentiel du sol:
1677	c
1678	DO k = 1, klev
1679	DO i = 1, klon
1680	zphi(i,k) = pphi(i,k) + pphis(i)
1681	ENDDO
1682	ENDDO
1683	c
1684	c Verifier les temperatures
1685	c
1686	cIM BEG
1687	IF (check) THEN
1688	amn=MIN(ftsol(1,is_ter),1000.)
1689	amx=MAX(ftsol(1,is_ter),-1000.)
1690	DO i=2, klon
1691	amn=MIN(ftsol(i,is_ter),amn)
1692	amx=MAX(ftsol(i,is_ter),amx)
1693	ENDDO
1694	c
1695	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1696	ENDIF !(check) THEN
1697	cIM END
1698	c
1699	CALL hgardfou(t_seri,ftsol,'debutphy')
1700	c
1701	cIM BEG
1702	IF (check) THEN
1703	amn=MIN(ftsol(1,is_ter),1000.)
1704	amx=MAX(ftsol(1,is_ter),-1000.)
1705	DO i=2, klon
1706	amn=MIN(ftsol(i,is_ter),amn)
1707	amx=MAX(ftsol(i,is_ter),amx)
1708	ENDDO
1709	c
1710	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1711	ENDIF !(check) THEN
1712	cIM END
1713	c
1714	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1715	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1716	c
1717	if (read_climoz >= 1) then
1718	C Ozone from a file
1719	! Update required ozone index:
1720	ro3i = int((days_elapsed + jh_cur - jh_1jan)
1721	$ / ioget_year_len(year_cur) * 360.) + 1
1722	if (ro3i == 361) ro3i = 360
1723	C (This should never occur, except perhaps because of roundup
1724	C error. See documentation.)
1725	if (ro3i /= co3i) then
1726	C Update ozone field:
1727	if (read_climoz == 1) then
1728	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1729	$ press_in_edg=press_climoz, paprs=paprs, v3=wo)
1730	else
1731	C read_climoz == 2
1732	call regr_pr_av(ncid_climoz,
1733	$ (/"tro3 ", "tro3_daylight"/),
1734	$ julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1735	$ v3=wo)
1736	end if
1737	! Convert from mole fraction of ozone to column density of ozone in a
1738	! cell, in kDU:
1739	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1740	$ * rmo3 / rmd * zmasse / dobson_u / 1e3
1741	C (By regridding ozone values for LMDZ only once every 360th of
1742	C year, we have already neglected the variation of pressure in one
1743	C 360th of year. So do not recompute "wo" at each time step even if
1744	C "zmasse" changes a little.)
1745	co3i = ro3i
1746	end if
1747	elseif (MOD(itap-1,lmt_pas) == 0) THEN
1748	C Once per day, update ozone from Royer:
1749	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1750	ENDIF
1751	c
1752	c Re-evaporer l'eau liquide nuageuse
1753	c
1754	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1755	DO i = 1, klon
1756	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1757	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1758	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1759	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1760	zb = MAX(0.0,ql_seri(i,k))
1761	za = - MAX(0.0,ql_seri(i,k))
1762	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1763	t_seri(i,k) = t_seri(i,k) + za
1764	q_seri(i,k) = q_seri(i,k) + zb
1765	ql_seri(i,k) = 0.0
1766	d_t_eva(i,k) = za
1767	d_q_eva(i,k) = zb
1768	ENDDO
1769	ENDDO
1770	cIM
1771	IF (ip_ebil_phy.ge.2) THEN
1772	ztit='after reevap'
1773	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1774	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1775	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1776	call diagphy(airephy,ztit,ip_ebil_phy
1777	e , zero_v, zero_v, zero_v, zero_v, zero_v
1778	e , zero_v, zero_v, zero_v, ztsol
1779	e , d_h_vcol, d_qt, d_ec
1780	s , fs_bound, fq_bound )
1781	C
1782	END IF
1783
1784	c
1785	c=========================================================================
1786	! Calculs de l'orbite.
1787	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1788	! doit donc etre placé avant radlwsw et pbl_surface
1789
1790	! calcul selon la routine utilisee pour les planetes
1791	if (new_orbit) then
1792	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1793	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1794	! day_since_equinox = (jD_cur) - jD_eq
1795	call solarlong(day_since_equinox, zlongi, dist)
1796	else
1797	! calcul selon la routine utilisee pour l'AR4
1798	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1799	! solarlong0
1800	if (solarlong0<-999.) then
1801	CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
1802	else
1803	zlongi=solarlong0 ! longitude solaire vraie
1804	dist=1. ! distance au soleil / moyenne
1805	endif
1806	endif
1807	if(prt_level.ge.1) &
1808	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1809
1810	! Avec ou sans cycle diurne
1811	IF (cycle_diurne) THEN
1812	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1813	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1814	ELSE
1815	CALL angle(zlongi, rlat, fract, rmu0)
1816	ENDIF
1817
1818	if (mydebug) then
1819	call writefield_phy('u_seri',u_seri,llm)
1820	call writefield_phy('v_seri',v_seri,llm)
1821	call writefield_phy('t_seri',t_seri,llm)
1822	call writefield_phy('q_seri',q_seri,llm)
1823	endif
1824
1825	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1826	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1827	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1828	c turbulent du couche limit.
1829	c
1830	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1831	c ecriture des fihiers hist_XXXX.nc, ces sont :
1832	c qsol, zq2m, s_pblh, s_lcl,
1833	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1834	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1835	c zxrugs, zu10m, zv10m, fder,
1836	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1837	c frugs, agesno, fsollw, fsolsw,
1838	c d_ts, fevap, fluxlat, t2m,
1839	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1840	c
1841	c Certains ne sont pas utiliser du tout :
1842	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1843	c
1844
1845	CALL pbl_surface(
1846	e dtime, date0, itap, days_elapsed+1,
1847	e debut, lafin,
1848	e rlon, rlat, rugoro, rmu0,
1849	e rain_fall, snow_fall, solsw, sollw,
1850	e t_seri, q_seri, u_seri, v_seri,
1851	e pplay, paprs, pctsrf,
1852	+ ftsol, falb1, falb2, u10m, v10m,
1853	s sollwdown, cdragh, cdragm, u1, v1,
1854	s albsol1, albsol2, sens, evap,
1855	s zxtsol, zxfluxlat, zt2m, qsat2m,
1856	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1857	s coefh, slab_wfbils,
1858	d qsol, zq2m, s_pblh, s_lcl,
1859	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1860	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1861	d zxrugs, zu10m, zv10m, fder,
1862	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1863	d frugs, agesno, fsollw, fsolsw,
1864	d d_ts, fevap, fluxlat, t2m,
1865	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1866	- dsens, devap, zxsnow,
1867	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1868
1869
1870	!-----------------------------------------------------------------------------------------
1871	! ajout des tendances de la diffusion turbulente
1872	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1873	!-----------------------------------------------------------------------------------------
1874
1875	if (mydebug) then
1876	call writefield_phy('u_seri',u_seri,llm)
1877	call writefield_phy('v_seri',v_seri,llm)
1878	call writefield_phy('t_seri',t_seri,llm)
1879	call writefield_phy('q_seri',q_seri,llm)
1880	endif
1881
1882
1883	IF (ip_ebil_phy.ge.2) THEN
1884	ztit='after surface_main'
1885	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1886	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1887	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1888	call diagphy(airephy,ztit,ip_ebil_phy
1889	e , zero_v, zero_v, zero_v, zero_v, sens
1890	e , evap , zero_v, zero_v, ztsol
1891	e , d_h_vcol, d_qt, d_ec
1892	s , fs_bound, fq_bound )
1893	END IF
1894
1895	c =================================================================== c
1896	c Calcul de Qsat
1897
1898	DO k = 1, klev
1899	DO i = 1, klon
1900	zx_t = t_seri(i,k)
1901	IF (thermcep) THEN
1902	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1903	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1904	zx_qs = MIN(0.5,zx_qs)
1905	zcor = 1./(1.-retv*zx_qs)
1906	zx_qs = zx_qs*zcor
1907	ELSE
1908	IF (zx_t.LT.t_coup) THEN
1909	zx_qs = qsats(zx_t)/pplay(i,k)
1910	ELSE
1911	zx_qs = qsatl(zx_t)/pplay(i,k)
1912	ENDIF
1913	ENDIF
1914	zqsat(i,k)=zx_qs
1915	ENDDO
1916	ENDDO
1917
1918	if (prt_level.ge.1) then
1919	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1920	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1921	endif
1922	c
1923	c Appeler la convection (au choix)
1924	c
1925	DO k = 1, klev
1926	DO i = 1, klon
1927	conv_q(i,k) = d_q_dyn(i,k)
1928	. + d_q_vdf(i,k)/dtime
1929	conv_t(i,k) = d_t_dyn(i,k)
1930	. + d_t_vdf(i,k)/dtime
1931	ENDDO
1932	ENDDO
1933	IF (check) THEN
1934	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1935	WRITE(lunout,*) "avantcon=", za
1936	ENDIF
1937	zx_ajustq = .FALSE.
1938	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1939	IF (zx_ajustq) THEN
1940	DO i = 1, klon
1941	z_avant(i) = 0.0
1942	ENDDO
1943	DO k = 1, klev
1944	DO i = 1, klon
1945	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1946	. *(paprs(i,k)-paprs(i,k+1))/RG
1947	ENDDO
1948	ENDDO
1949	ENDIF
1950
1951	c Calcule de vitesse verticale a partir de flux de masse verticale
1952	DO k = 1, klev
1953	DO i = 1, klon
1954	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1955	END DO
1956	END DO
1957	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
1958	$ omega(igout, :)
1959
1960	IF (iflag_con.EQ.1) THEN
1961	stop'reactiver le call conlmd dans physiq.F'
1962	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1963	c . d_t_con, d_q_con,
1964	c . rain_con, snow_con, ibas_con, itop_con)
1965	ELSE IF (iflag_con.EQ.2) THEN
1966	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1967	e conv_t, conv_q, -evap, omega,
1968	s d_t_con, d_q_con, rain_con, snow_con,
1969	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1970	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1971	d_u_con = 0.
1972	d_v_con = 0.
1973
1974	WHERE (rain_con < 0.) rain_con = 0.
1975	WHERE (snow_con < 0.) snow_con = 0.
1976	DO i = 1, klon
1977	ibas_con(i) = klev+1 - kcbot(i)
1978	itop_con(i) = klev+1 - kctop(i)
1979	ENDDO
1980	ELSE IF (iflag_con.GE.3) THEN
1981	c nb of tracers for the KE convection:
1982	c MAF la partie traceurs est faite dans phytrac
1983	c on met ntra=1 pour limiter les appels mais on peut
1984	c supprimer les calculs / ftra.
1985	ntra = 1
1986
1987	c=====================================================================================
1988	cajout pour la parametrisation des poches froides:
1989	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1990	do k=1,klev
1991	do i=1,klon
1992	if (iflag_wake.eq.1) then
1993	t_wake(i,k) = t_seri(i,k)
1994	. +(1-wake_s(i))*wake_deltat(i,k)
1995	q_wake(i,k) = q_seri(i,k)
1996	. +(1-wake_s(i))*wake_deltaq(i,k)
1997	t_undi(i,k) = t_seri(i,k)
1998	. -wake_s(i)*wake_deltat(i,k)
1999	q_undi(i,k) = q_seri(i,k)
2000	. -wake_s(i)*wake_deltaq(i,k)
2001	else
2002	t_wake(i,k) = t_seri(i,k)
2003	q_wake(i,k) = q_seri(i,k)
2004	t_undi(i,k) = t_seri(i,k)
2005	q_undi(i,k) = q_seri(i,k)
2006	endif
2007	enddo
2008	enddo
2009
2010	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2011	cc-- pour le soulevement des particules dans le modele convectif
2012	c
2013	do i = 1,klon
2014	ALE(i) = 0.
2015	ALP(i) = 0.
2016	enddo
2017	c
2018	ccalcul de ale_wake et alp_wake
2019	do i = 1,klon
2020	if (iflag_wake.eq.1) then
2021	ale_wake(i) = 0.5wake_cstar(i)*2
2022	alp_wake(i) = wake_fip(i)
2023	else
2024	ale_wake(i) = 0.
2025	alp_wake(i) = 0.
2026	endif
2027	enddo
2028	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2029	cdans le thermique sinon
2030	if (iflag_coupl.eq.0) then
2031	if (debut) print*,'ALE et ALP imposes'
2032	do i = 1,klon
2033	con ne couple que ale
2034	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2035	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2036	con ne couple que alp
2037	c ALP(i) = alp_wake(i) + Alp_bl(i)
2038	ALP(i) = alp_wake(i) + alp_bl_prescr
2039	enddo
2040	else
2041	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2042	do i = 1,klon
2043	ALE(i) = max(ale_wake(i),Ale_bl(i))
2044	ALP(i) = alp_wake(i) + Alp_bl(i)
2045	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2046	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2047	enddo
2048	endif
2049	do i=1,klon
2050	if (alp(i)>alp_max) then
2051	IF(prt_level>9)WRITE(lunout,*) &
2052	& 'WARNING SUPER ALP (seuil=',alp_max,
2053	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2054	alp(i)=alp_max
2055	endif
2056	if (ale(i)>ale_max) then
2057	IF(prt_level>9)WRITE(lunout,*) &
2058	& 'WARNING SUPER ALE (seuil=',ale_max,
2059	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2060	ale(i)=ale_max
2061	endif
2062	enddo
2063
2064	cfin calcul ale et alp
2065	c=================================================================================================
2066
2067
2068	c sb, oct02:
2069	c Schema de convection modularise et vectorise:
2070	c (driver commun aux versions 3 et 4)
2071	c
2072	IF (ok_cvl) THEN ! new driver for convectL
2073
2074	CALL concvl (iflag_con,iflag_clos,
2075	. dtime,paprs,pplay,t_undi,q_undi,
2076	. t_wake,q_wake,wake_s,
2077	. u_seri,v_seri,tr_seri,nbtr,
2078	. ALE,ALP,
2079	. ema_work1,ema_work2,
2080	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2081	. rain_con, snow_con, ibas_con, itop_con, sigd,
2082	. upwd,dnwd,dnwd0,
2083	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2084	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2085	. pmflxr,pmflxs,da,phi,mp,
2086	. ftd,fqd,lalim_conv,wght_th)
2087
2088	cIM begin
2089	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2090	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2091	cIM end
2092	cIM cf. FH
2093	clwcon0=qcondc
2094	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2095
2096	ELSE ! ok_cvl
2097	c MAF conema3 ne contient pas les traceurs
2098	CALL conema3 (dtime,
2099	. paprs,pplay,t_seri,q_seri,
2100	. u_seri,v_seri,tr_seri,ntra,
2101	. ema_work1,ema_work2,
2102	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2103	. rain_con, snow_con, ibas_con, itop_con,
2104	. upwd,dnwd,dnwd0,bas,top,
2105	. Ma,cape,tvp,rflag,
2106	. pbase
2107	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2108	. ,clwcon0)
2109
2110	ENDIF ! ok_cvl
2111
2112	c
2113	c Correction precip
2114	rain_con = rain_con * cvl_corr
2115	snow_con = snow_con * cvl_corr
2116	c
2117
2118	IF (.NOT. ok_gust) THEN
2119	do i = 1, klon
2120	wd(i)=0.0
2121	enddo
2122	ENDIF
2123
2124	c =================================================================== c
2125	c Calcul des proprietes des nuages convectifs
2126	c
2127
2128	c calcul des proprietes des nuages convectifs
2129	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2130	call clouds_gno
2131	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2132
2133	c =================================================================== c
2134
2135	DO i = 1, klon
2136	ema_pcb(i) = pbase(i)
2137	ENDDO
2138	DO i = 1, klon
2139
2140	! L'idicage de itop_con peut cacher un pb potentiel
2141	! FH sous la dictee de JYG, CR
2142	ema_pct(i) = paprs(i,itop_con(i)+1)
2143
2144	if (itop_con(i).gt.klev-3) then
2145	if(prt_level >= 9) then
2146	write(lunout,*)'La convection monte trop haut '
2147	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2148	endif
2149	endif
2150	ENDDO
2151	DO i = 1, klon
2152	ema_cbmf(i) = ema_workcbmf(i)
2153	ENDDO
2154	ELSE IF (iflag_con.eq.0) THEN
2155	write(lunout,*) 'On n appelle pas la convection'
2156	clwcon0=0.
2157	rnebcon0=0.
2158	d_t_con=0.
2159	d_q_con=0.
2160	d_u_con=0.
2161	d_v_con=0.
2162	rain_con=0.
2163	snow_con=0.
2164	bas=1
2165	top=1
2166	ELSE
2167	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2168	CALL abort
2169	ENDIF
2170
2171	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2172	c . d_u_con, d_v_con)
2173
2174	!-----------------------------------------------------------------------------------------
2175	! ajout des tendances de la diffusion turbulente
2176	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2177	!-----------------------------------------------------------------------------------------
2178
2179	if (mydebug) then
2180	call writefield_phy('u_seri',u_seri,llm)
2181	call writefield_phy('v_seri',v_seri,llm)
2182	call writefield_phy('t_seri',t_seri,llm)
2183	call writefield_phy('q_seri',q_seri,llm)
2184	endif
2185
2186	cIM
2187	IF (ip_ebil_phy.ge.2) THEN
2188	ztit='after convect'
2189	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2190	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2191	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2192	call diagphy(airephy,ztit,ip_ebil_phy
2193	e , zero_v, zero_v, zero_v, zero_v, zero_v
2194	e , zero_v, rain_con, snow_con, ztsol
2195	e , d_h_vcol, d_qt, d_ec
2196	s , fs_bound, fq_bound )
2197	END IF
2198	C
2199	IF (check) THEN
2200	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2201	WRITE(lunout,*)"aprescon=", za
2202	zx_t = 0.0
2203	za = 0.0
2204	DO i = 1, klon
2205	za = za + airephy(i)/FLOAT(klon)
2206	zx_t = zx_t + (rain_con(i)+
2207	. snow_con(i))*airephy(i)/FLOAT(klon)
2208	ENDDO
2209	zx_t = zx_t/za*dtime
2210	WRITE(lunout,*)"Precip=", zx_t
2211	ENDIF
2212	IF (zx_ajustq) THEN
2213	DO i = 1, klon
2214	z_apres(i) = 0.0
2215	ENDDO
2216	DO k = 1, klev
2217	DO i = 1, klon
2218	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2219	. *(paprs(i,k)-paprs(i,k+1))/RG
2220	ENDDO
2221	ENDDO
2222	DO i = 1, klon
2223	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2224	. /z_apres(i)
2225	ENDDO
2226	DO k = 1, klev
2227	DO i = 1, klon
2228	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2229	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2230	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2231	ENDIF
2232	ENDDO
2233	ENDDO
2234	ENDIF
2235	zx_ajustq=.FALSE.
2236
2237	c
2238	c=============================================================================
2239	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2240	cpour la couche limite diffuse pour l instant
2241	c
2242	if (iflag_wake.eq.1) then
2243	DO k=1,klev
2244	DO i=1,klon
2245	dt_dwn(i,k) = ftd(i,k)
2246	wdt_PBL(i,k) = 0.
2247	dq_dwn(i,k) = fqd(i,k)
2248	wdq_PBL(i,k) = 0.
2249	M_dwn(i,k) = dnwd0(i,k)
2250	M_up(i,k) = upwd(i,k)
2251	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2252	udt_PBL(i,k) = 0.
2253	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2254	udq_PBL(i,k) = 0.
2255	ENDDO
2256	ENDDO
2257	c
2258	ccalcul caracteristiques de la poche froide
2259	call calWAKE (paprs,pplay,dtime
2260	: ,t_seri,q_seri,omega
2261	: ,dt_dwn,dq_dwn,M_dwn,M_up
2262	: ,dt_a,dq_a,sigd
2263	: ,wdt_PBL,wdq_PBL
2264	: ,udt_PBL,udq_PBL
2265	o ,wake_deltat,wake_deltaq,wake_dth
2266	o ,wake_h,wake_s,wake_dens
2267	o ,wake_pe,wake_fip,wake_gfl
2268	o ,dt_wake,dq_wake
2269	o ,wake_k, t_undi,q_undi
2270	o ,wake_omgbdth,wake_dp_omgb
2271	o ,wake_dtKE,wake_dqKE
2272	o ,wake_dtPBL,wake_dqPBL
2273	o ,wake_omg,wake_dp_deltomg
2274	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2275	o ,wake_ddeltat,wake_ddeltaq)
2276	c
2277	!-----------------------------------------------------------------------------------------
2278	! ajout des tendances des poches froides
2279	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2280	! coherent avec les autres d_t_...
2281	d_t_wake(:,:)=dt_wake(:,:)*dtime
2282	d_q_wake(:,:)=dq_wake(:,:)*dtime
2283	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2284	!-----------------------------------------------------------------------------------------
2285
2286	endif
2287	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2288	c
2289	c===================================================================
2290	c Convection seche (thermiques ou ajustement)
2291	c===================================================================
2292	c
2293	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2294	s ,seuil_inversion,weak_inversion,dthmin)
2295
2296
2297
2298	d_t_ajsb(:,:)=0.
2299	d_q_ajsb(:,:)=0.
2300	d_t_ajs(:,:)=0.
2301	d_u_ajs(:,:)=0.
2302	d_v_ajs(:,:)=0.
2303	d_q_ajs(:,:)=0.
2304	clwcon0th(:,:)=0.
2305	c
2306	fm_therm(:,:)=0.
2307	entr_therm(:,:)=0.
2308	detr_therm(:,:)=0.
2309	c
2310	IF(prt_level>9)WRITE(lunout,*)
2311	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2312	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2313	if(iflag_thermals.lt.0) then
2314	c Rien
2315	c ====
2316	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2317
2318
2319	else
2320
2321	c Thermiques
2322	c ==========
2323	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2324	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2325
2326
2327	if (iflag_thermals.gt.1) then
2328	call calltherm(pdtphys
2329	s ,pplay,paprs,pphi,weak_inversion
2330	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2331	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2332	s ,fm_therm,entr_therm,detr_therm
2333	s ,zqasc,clwcon0th,lmax_th,ratqscth
2334	s ,ratqsdiff,zqsatth
2335	con rajoute ale et alp, et les caracteristiques de la couche alim
2336	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
2337	endif
2338
2339
2340	c Ajustement sec
2341	c ==============
2342
2343	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2344	! a partir du sommet des thermiques.
2345	! Dans le cas contraire, on demarre au niveau 1.
2346
2347	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2348
2349	if(iflag_thermals.eq.0) then
2350	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2351	limbas(:)=1
2352	else
2353	limbas(:)=lmax_th(:)
2354	endif
2355
2356	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2357	! pour des test de convergence numerique.
2358	! Le nouveau ajsec est a priori mieux, meme pour le cas
2359	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2360	! non nulles numeriquement pour des mailles non concernees.
2361
2362	if (iflag_thermals.eq.0) then
2363	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2364	s , d_t_ajsb, d_q_ajsb)
2365	else
2366	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2367	s , d_t_ajsb, d_q_ajsb)
2368	endif
2369
2370	!-----------------------------------------------------------------------------------------
2371	! ajout des tendances de l'ajustement sec ou des thermiques
2372	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2373	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2374	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2375
2376	!-----------------------------------------------------------------------------------------
2377
2378	endif
2379
2380	endif
2381	c
2382	c===================================================================
2383	cIM
2384	IF (ip_ebil_phy.ge.2) THEN
2385	ztit='after dry_adjust'
2386	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2387	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2388	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2389	END IF
2390
2391
2392	c-------------------------------------------------------------------------
2393	c Caclul des ratqs
2394	c-------------------------------------------------------------------------
2395
2396	c print*,'calcul des ratqs'
2397	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2398	c ----------------
2399	c on ecrase le tableau ratqsc calcule par clouds_gno
2400	if (iflag_cldcon.eq.1) then
2401	do k=1,klev
2402	do i=1,klon
2403	if(ptconv(i,k)) then
2404	ratqsc(i,k)=ratqsbas
2405	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2406	else
2407	ratqsc(i,k)=0.
2408	endif
2409	enddo
2410	enddo
2411
2412	c-----------------------------------------------------------------------
2413	c par nversion de la fonction log normale
2414	c-----------------------------------------------------------------------
2415	else if (iflag_cldcon.eq.4) then
2416	ptconvth(:,:)=.false.
2417	ratqsc(:,:)=0.
2418	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2419	call clouds_gno
2420	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2421	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2422
2423	c-----------------------------------------------------------------------
2424	c par calcul direct de l'ecart-type
2425	c-----------------------------------------------------------------------
2426
2427	else if (iflag_cldcon>=5) then
2428	wmax_th(:)=0.
2429	zmax_th(:)=0.
2430	do k=1,klev
2431	do i=1,klon
2432	wmax_th(i)=max(wmax_th(i),zw2(i,k))
2433	if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
2434	enddo
2435	enddo
2436	tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
2437	print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
2438
2439	c On impose que l'air autour de la fraction couverte par le thermique
2440	c plus son air detraine durant tau_overturning_th soit superieur
2441	c a 0.1 q_seri
2442	zz=0.1
2443	do k=1,klev
2444	do i=1,klon
2445	lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
2446	s tau_overturning_th(i)*detr_therm(i,k)
2447	s *rg/(paprs(i,k)-paprs(i,k+1))
2448	znum=(1.-zz)*q_seri(i,k)
2449	zden=zqasc(i,k)-zz*q_seri(i,k)
2450	if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
2451	lambda_th(i,k)=min(lambda_th(i,k),0.9)
2452	enddo
2453	enddo
2454
2455	if(iflag_cldcon==5) then
2456	do k=1,klev
2457	do i=1,klon
2458	ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
2459	s abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
2460	enddo
2461	enddo
2462	else if(iflag_cldcon==6) then
2463	do k=1,klev
2464	do i=1,klon
2465	ratqsc(i,k)=sqrt(lambda_th(i,k))*
2466	s (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
2467	enddo
2468	enddo
2469	endif
2470
2471	endif
2472
2473	c ratqs stables
2474	c -------------
2475
2476	if (iflag_ratqs.eq.0) then
2477
2478	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2479	do k=1,klev
2480	do i=1, klon
2481	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2482	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2483	enddo
2484	enddo
2485
2486	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2487	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2488	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2489	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2490	! Il s'agit de differents tests dans la phase de reglage du modele
2491	! avec thermiques.
2492
2493	else if (iflag_ratqs.eq.1) then
2494
2495	do k=1,klev
2496	do i=1, klon
2497	if (pplay(i,k).ge.60000.) then
2498	ratqss(i,k)=ratqsbas
2499	else if ((pplay(i,k).ge.30000.).and.
2500	s (pplay(i,k).lt.60000.)) then
2501	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2502	s (60000.-pplay(i,k))/(60000.-30000.)
2503	else
2504	ratqss(i,k)=ratqshaut
2505	endif
2506	enddo
2507	enddo
2508
2509	else
2510
2511	do k=1,klev
2512	do i=1, klon
2513	if (pplay(i,k).ge.60000.) then
2514	ratqss(i,k)=ratqsbas
2515	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2516	else if ((pplay(i,k).ge.30000.).and.
2517	s (pplay(i,k).lt.60000.)) then
2518	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2519	s (60000.-pplay(i,k))/(60000.-30000.)
2520	else
2521	ratqss(i,k)=ratqshaut
2522	endif
2523	enddo
2524	enddo
2525	endif
2526
2527
2528
2529
2530	c ratqs final
2531	c -----------
2532
2533	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2534	s .or.iflag_cldcon.ge.4) then
2535
2536	! On ajoute une constante au ratqsc*2 pour tenir compte de
2537	! fluctuations turbulentes de petite echelle
2538
2539	do k=1,klev
2540	do i=1,klon
2541	if ((fm_therm(i,k).gt.1.e-10)) then
2542	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2543	endif
2544	enddo
2545	enddo
2546
2547	! les ratqs sont une combinaison de ratqss et ratqsc
2548	print*,'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
2549
2550	if (tau_ratqs>1.e-10) then
2551	facteur=exp(-pdtphys/tau_ratqs)
2552	else
2553	facteur=0.
2554	endif
2555	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2556	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2557	! FH 22/09/2009
2558	! La ligne ci-dessous faisait osciller le modele et donnait une solution
2559	! assymptotique bidon et dÃ©pendant fortement du pas de temps.
2560	! ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2561	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2562	ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
2563	else
2564	! on ne prend que le ratqs stable pour fisrtilp
2565	ratqs(:,:)=ratqss(:,:)
2566	endif
2567
2568
2569	c
2570	c Appeler le processus de condensation a grande echelle
2571	c et le processus de precipitation
2572	c-------------------------------------------------------------------------
2573	CALL fisrtilp(dtime,paprs,pplay,
2574	. t_seri, q_seri,ptconv,ratqs,
2575	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2576	. rain_lsc, snow_lsc,
2577	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2578	. frac_impa, frac_nucl,
2579	. prfl, psfl, rhcl)
2580
2581	WHERE (rain_lsc < 0) rain_lsc = 0.
2582	WHERE (snow_lsc < 0) snow_lsc = 0.
2583	!-----------------------------------------------------------------------------------------
2584	! ajout des tendances de la diffusion turbulente
2585	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2586	!-----------------------------------------------------------------------------------------
2587	DO k = 1, klev
2588	DO i = 1, klon
2589	cldfra(i,k) = rneb(i,k)
2590	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2591	ENDDO
2592	ENDDO
2593	IF (check) THEN
2594	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2595	WRITE(lunout,*)"apresilp=", za
2596	zx_t = 0.0
2597	za = 0.0
2598	DO i = 1, klon
2599	za = za + airephy(i)/FLOAT(klon)
2600	zx_t = zx_t + (rain_lsc(i)
2601	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2602	ENDDO
2603	zx_t = zx_t/za*dtime
2604	WRITE(lunout,*)"Precip=", zx_t
2605	ENDIF
2606	cIM
2607	IF (ip_ebil_phy.ge.2) THEN
2608	ztit='after fisrt'
2609	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2610	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2611	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2612	call diagphy(airephy,ztit,ip_ebil_phy
2613	e , zero_v, zero_v, zero_v, zero_v, zero_v
2614	e , zero_v, rain_lsc, snow_lsc, ztsol
2615	e , d_h_vcol, d_qt, d_ec
2616	s , fs_bound, fq_bound )
2617	END IF
2618
2619	if (mydebug) then
2620	call writefield_phy('u_seri',u_seri,llm)
2621	call writefield_phy('v_seri',v_seri,llm)
2622	call writefield_phy('t_seri',t_seri,llm)
2623	call writefield_phy('q_seri',q_seri,llm)
2624	endif
2625
2626	c
2627	c-------------------------------------------------------------------
2628	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2629	c-------------------------------------------------------------------
2630
2631	c 1. NUAGES CONVECTIFS
2632	c
2633	cIM cf FH
2634	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2635	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2636	snow_tiedtke=0.
2637	c print*,'avant calcul de la pseudo precip '
2638	c print*,'iflag_cldcon',iflag_cldcon
2639	if (iflag_cldcon.eq.-1) then
2640	rain_tiedtke=rain_con
2641	else
2642	c print*,'calcul de la pseudo precip '
2643	rain_tiedtke=0.
2644	c print*,'calcul de la pseudo precip 0'
2645	do k=1,klev
2646	do i=1,klon
2647	if (d_q_con(i,k).lt.0.) then
2648	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2649	s *(paprs(i,k)-paprs(i,k+1))/rg
2650	endif
2651	enddo
2652	enddo
2653	endif
2654	c
2655	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2656	c
2657
2658	c Nuages diagnostiques pour Tiedtke
2659	CALL diagcld1(paprs,pplay,
2660	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2661	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2662	. diafra,dialiq)
2663	DO k = 1, klev
2664	DO i = 1, klon
2665	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2666	cldliq(i,k) = dialiq(i,k)
2667	cldfra(i,k) = diafra(i,k)
2668	ENDIF
2669	ENDDO
2670	ENDDO
2671
2672	ELSE IF (iflag_cldcon.ge.3) THEN
2673	c On prend pour les nuages convectifs le max du calcul de la
2674	c convection et du calcul du pas de temps precedent diminue d'un facteur
2675	c facttemps
2676	facteur = pdtphys *facttemps
2677	do k=1,klev
2678	do i=1,klon
2679	rnebcon(i,k)=rnebcon(i,k)*facteur
2680	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2681	s then
2682	rnebcon(i,k)=rnebcon0(i,k)
2683	clwcon(i,k)=clwcon0(i,k)
2684	endif
2685	enddo
2686	enddo
2687
2688	c
2689	cjq - introduce the aerosol direct and first indirect radiative forcings
2690	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2691	IF (ok_ade.OR.ok_aie) THEN
2692	IF (.NOT. aerosol_couple)
2693	& CALL readaerosol_optic(
2694	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
2695	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
2696	& mass_solu_aero, mass_solu_aero_pi,
2697	& tau_aero, piz_aero, cg_aero,
2698	& tausum_aero, tau3d_aero)
2699	ELSE
2700	tau_aero(:,:,:,:) = 0.
2701	piz_aero(:,:,:,:) = 0.
2702	cg_aero(:,:,:,:) = 0.
2703	ENDIF
2704
2705	cIM calcul nuages par le simulateur ISCCP
2706	c
2707	#ifdef histISCCP
2708	IF (ok_isccp) THEN
2709	c
2710	cIM lecture invtau, tautab des fichiers formattes
2711	c
2712	IF (debut) THEN
2713	c$OMP MASTER
2714	c
2715	open(99,file='tautab.formatted', FORM='FORMATTED')
2716	read(99,'(f30.20)') tautab_omp
2717	close(99)
2718	c
2719	open(99,file='invtau.formatted',form='FORMATTED')
2720	read(99,'(i10)') invtau_omp
2721
2722	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
2723	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2724
2725	close(99)
2726	c$OMP END MASTER
2727	c$OMP BARRIER
2728	tautab=tautab_omp
2729	invtau=invtau_omp
2730	c
2731	ENDIF !debut
2732	c
2733	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2734	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2735	#include "calcul_simulISCCP.h"
2736	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2737	ENDIF !ok_isccp
2738	#endif
2739
2740	c On prend la somme des fractions nuageuses et des contenus en eau
2741	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2742	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2743
2744	ENDIF
2745	c
2746	c 2. NUAGES STARTIFORMES
2747	c
2748	IF (ok_stratus) THEN
2749	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2750	DO k = 1, klev
2751	DO i = 1, klon
2752	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2753	cldliq(i,k) = dialiq(i,k)
2754	cldfra(i,k) = diafra(i,k)
2755	ENDIF
2756	ENDDO
2757	ENDDO
2758	ENDIF
2759	c
2760	c Precipitation totale
2761	c
2762	DO i = 1, klon
2763	rain_fall(i) = rain_con(i) + rain_lsc(i)
2764	snow_fall(i) = snow_con(i) + snow_lsc(i)
2765	ENDDO
2766	cIM
2767	IF (ip_ebil_phy.ge.2) THEN
2768	ztit="after diagcld"
2769	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2770	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2771	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2772	END IF
2773	c
2774	c Calculer l'humidite relative pour diagnostique
2775	c
2776	DO k = 1, klev
2777	DO i = 1, klon
2778	zx_t = t_seri(i,k)
2779	IF (thermcep) THEN
2780	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2781	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2782	zx_qs = MIN(0.5,zx_qs)
2783	zcor = 1./(1.-retv*zx_qs)
2784	zx_qs = zx_qs*zcor
2785	ELSE
2786	IF (zx_t.LT.t_coup) THEN
2787	zx_qs = qsats(zx_t)/pplay(i,k)
2788	ELSE
2789	zx_qs = qsatl(zx_t)/pplay(i,k)
2790	ENDIF
2791	ENDIF
2792	zx_rh(i,k) = q_seri(i,k)/zx_qs
2793	zqsat(i,k)=zx_qs
2794	ENDDO
2795	ENDDO
2796
2797	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2798	c equivalente a 2m (tpote) pour diagnostique
2799	c
2800	DO i = 1, klon
2801	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2802	IF (thermcep) THEN
2803	IF(zt2m(i).LT.RTT) then
2804	Lheat=RLSTT
2805	ELSE
2806	Lheat=RLVTT
2807	ENDIF
2808	ELSE
2809	IF (zt2m(i).LT.RTT) THEN
2810	Lheat=RLSTT
2811	ELSE
2812	Lheat=RLVTT
2813	ENDIF
2814	ENDIF
2815	tpote(i) = tpot(i)*
2816	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2817	ENDDO
2818
2819	IF (config_inca /= 'none') THEN
2820	#ifdef INCA
2821	CALL VTe(VTphysiq)
2822	CALL VTb(VTinca)
2823	calday = FLOAT(days_elapsed + 1) + jH_cur
2824
2825	IF (config_inca == 'aero') THEN
2826	CALL AEROSOL_METEO_CALC(
2827	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2828	$ prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
2829	END IF
2830
2831	zxsnow_dummy(:) = 0.0
2832
2833	CALL chemhook_begin (calday,
2834	$ days_elapsed+1,
2835	$ jH_cur,
2836	$ pctsrf(1,1),
2837	$ rlat,
2838	$ rlon,
2839	$ airephy,
2840	$ paprs,
2841	$ pplay,
2842	$ coefh,
2843	$ pphi,
2844	$ t_seri,
2845	$ u,
2846	$ v,
2847	$ wo(:, :, 1),
2848	$ q_seri,
2849	$ zxtsol,
2850	$ zxsnow_dummy,
2851	$ solsw,
2852	$ albsol1,
2853	$ rain_fall,
2854	$ snow_fall,
2855	$ itop_con,
2856	$ ibas_con,
2857	$ cldfra,
2858	$ iim,
2859	$ jjm,
2860	$ tr_seri,
2861	$ ftsol,
2862	$ paprs,
2863	$ cdragh,
2864	$ cdragm,
2865	$ pctsrf,
2866	$ pdtphys,
2867	$ itap)
2868
2869	CALL VTe(VTinca)
2870	CALL VTb(VTphysiq)
2871	#endif
2872	END IF !config_inca /= 'none'
2873	c
2874	c Calculer les parametres optiques des nuages et quelques
2875	c parametres pour diagnostiques:
2876	c
2877
2878	IF (aerosol_couple) THEN
2879	mass_solu_aero(:,:) = ccm(:,:,1)
2880	mass_solu_aero_pi(:,:) = ccm(:,:,2)
2881	END IF
2882
2883	if (ok_newmicro) then
2884	CALL newmicro (paprs, pplay,ok_newmicro,
2885	. t_seri, cldliq, cldfra, cldtau, cldemi,
2886	. cldh, cldl, cldm, cldt, cldq,
2887	. flwp, fiwp, flwc, fiwc,
2888	e ok_aie,
2889	e mass_solu_aero, mass_solu_aero_pi,
2890	e bl95_b0, bl95_b1,
2891	s cldtaupi, re, fl, ref_liq, ref_ice)
2892	else
2893	CALL nuage (paprs, pplay,
2894	. t_seri, cldliq, cldfra, cldtau, cldemi,
2895	. cldh, cldl, cldm, cldt, cldq,
2896	e ok_aie,
2897	e mass_solu_aero, mass_solu_aero_pi,
2898	e bl95_b0, bl95_b1,
2899	s cldtaupi, re, fl)
2900
2901	endif
2902	c
2903	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2904	c
2905	IF (MOD(itaprad,radpas).EQ.0) THEN
2906
2907	DO i = 1, klon
2908	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2909	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2910	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2911	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2912	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2913	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2914	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2915	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2916	ENDDO
2917
2918	if (mydebug) then
2919	call writefield_phy('u_seri',u_seri,llm)
2920	call writefield_phy('v_seri',v_seri,llm)
2921	call writefield_phy('t_seri',t_seri,llm)
2922	call writefield_phy('q_seri',q_seri,llm)
2923	endif
2924
2925	IF (aerosol_couple) THEN
2926	#ifdef INCA
2927	CALL radlwsw_inca
2928	e (kdlon,kflev,dist, rmu0, fract, solaire,
2929	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2930	e wo(:, :, 1),
2931	e cldfra, cldemi, cldtau,
2932	s heat,heat0,cool,cool0,radsol,albpla,
2933	s topsw,toplw,solsw,sollw,
2934	s sollwdown,
2935	s topsw0,toplw0,solsw0,sollw0,
2936	s lwdn0, lwdn, lwup0, lwup,
2937	s swdn0, swdn, swup0, swup,
2938	e ok_ade, ok_aie,
2939	e tau_aero, piz_aero, cg_aero,
2940	s topswad_aero, solswad_aero,
2941	s topswad0_aero, solswad0_aero,
2942	s topsw_aero, topsw0_aero,
2943	s solsw_aero, solsw0_aero,
2944	e cldtaupi,
2945	s topswai_aero, solswai_aero)
2946
2947	#endif
2948	ELSE
2949
2950	CALL radlwsw
2951	e (dist, rmu0, fract,
2952	e paprs, pplay,zxtsol,albsol1, albsol2,
2953	e t_seri,q_seri,wo,
2954	e cldfra, cldemi, cldtau,
2955	e ok_ade, ok_aie,
2956	e tau_aero, piz_aero, cg_aero,
2957	e cldtaupi,new_aod,
2958	e zqsat, flwc, fiwc,
2959	s heat,heat0,cool,cool0,radsol,albpla,
2960	s topsw,toplw,solsw,sollw,
2961	s sollwdown,
2962	s topsw0,toplw0,solsw0,sollw0,
2963	s lwdn0, lwdn, lwup0, lwup,
2964	s swdn0, swdn, swup0, swup,
2965	s topswad_aero, solswad_aero,
2966	s topswai_aero, solswai_aero,
2967	o topswad0_aero, solswad0_aero,
2968	o topsw_aero, topsw0_aero,
2969	o solsw_aero, solsw0_aero,
2970	o topswcf_aero, solswcf_aero)
2971
2972
2973	ENDIF ! aerosol_couple
2974	itaprad = 0
2975	ENDIF ! MOD(itaprad,radpas)
2976	itaprad = itaprad + 1
2977
2978	if (iflag_radia.eq.0) then
2979	print *,'--------------------------------------------------'
2980	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
2981	print *,'>>>> heat et cool mis a zero '
2982	print *,'--------------------------------------------------'
2983	heat=0.
2984	cool=0.
2985	endif
2986
2987	c
2988	c Ajouter la tendance des rayonnements (tous les pas)
2989	c
2990	DO k = 1, klev
2991	DO i = 1, klon
2992	t_seri(i,k) = t_seri(i,k)
2993	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
2994	ENDDO
2995	ENDDO
2996	c
2997	if (mydebug) then
2998	call writefield_phy('u_seri',u_seri,llm)
2999	call writefield_phy('v_seri',v_seri,llm)
3000	call writefield_phy('t_seri',t_seri,llm)
3001	call writefield_phy('q_seri',q_seri,llm)
3002	endif
3003
3004	cIM
3005	IF (ip_ebil_phy.ge.2) THEN
3006	ztit='after rad'
3007	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3008	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3009	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3010	call diagphy(airephy,ztit,ip_ebil_phy
3011	e , topsw, toplw, solsw, sollw, zero_v
3012	e , zero_v, zero_v, zero_v, ztsol
3013	e , d_h_vcol, d_qt, d_ec
3014	s , fs_bound, fq_bound )
3015	END IF
3016	c
3017	c
3018	c Calculer l'hydrologie de la surface
3019	c
3020	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3021	c . agesno, ftsol,fqsurf,fsnow, ruis)
3022	c
3023
3024	c
3025	c Calculer le bilan du sol et la derive de temperature (couplage)
3026	c
3027	DO i = 1, klon
3028	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3029	c a la demande de JLD
3030	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3031	ENDDO
3032	c
3033	cmoddeblott(jan95)
3034	c Appeler le programme de parametrisation de l'orographie
3035	c a l'echelle sous-maille:
3036	c
3037	IF (ok_orodr) THEN
3038	c
3039	c selection des points pour lesquels le shema est actif:
3040	igwd=0
3041	DO i=1,klon
3042	itest(i)=0
3043	c IF ((zstd(i).gt.10.0)) THEN
3044	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3045	itest(i)=1
3046	igwd=igwd+1
3047	idx(igwd)=i
3048	ENDIF
3049	ENDDO
3050	c igwdim=MAX(1,igwd)
3051	c
3052	IF (ok_strato) THEN
3053
3054	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3055	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3056	e igwd,idx,itest,
3057	e t_seri, u_seri, v_seri,
3058	s zulow, zvlow, zustrdr, zvstrdr,
3059	s d_t_oro, d_u_oro, d_v_oro)
3060
3061	ELSE
3062	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3063	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3064	e igwd,idx,itest,
3065	e t_seri, u_seri, v_seri,
3066	s zulow, zvlow, zustrdr, zvstrdr,
3067	s d_t_oro, d_u_oro, d_v_oro)
3068	ENDIF
3069	c
3070	c ajout des tendances
3071	!-----------------------------------------------------------------------------------------
3072	! ajout des tendances de la trainee de l'orographie
3073	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3074	!-----------------------------------------------------------------------------------------
3075	c
3076	ENDIF ! fin de test sur ok_orodr
3077	c
3078	if (mydebug) then
3079	call writefield_phy('u_seri',u_seri,llm)
3080	call writefield_phy('v_seri',v_seri,llm)
3081	call writefield_phy('t_seri',t_seri,llm)
3082	call writefield_phy('q_seri',q_seri,llm)
3083	endif
3084
3085	IF (ok_orolf) THEN
3086	c
3087	c selection des points pour lesquels le shema est actif:
3088	igwd=0
3089	DO i=1,klon
3090	itest(i)=0
3091	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3092	itest(i)=1
3093	igwd=igwd+1
3094	idx(igwd)=i
3095	ENDIF
3096	ENDDO
3097	c igwdim=MAX(1,igwd)
3098	c
3099	IF (ok_strato) THEN
3100
3101	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3102	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3103	e igwd,idx,itest,
3104	e t_seri, u_seri, v_seri,
3105	s zulow, zvlow, zustrli, zvstrli,
3106	s d_t_lif, d_u_lif, d_v_lif )
3107
3108	ELSE
3109	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3110	e rlat,zmea,zstd,zpic,
3111	e itest,
3112	e t_seri, u_seri, v_seri,
3113	s zulow, zvlow, zustrli, zvstrli,
3114	s d_t_lif, d_u_lif, d_v_lif)
3115	ENDIF
3116	c
3117	!-----------------------------------------------------------------------------------------
3118	! ajout des tendances de la portance de l'orographie
3119	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3120	!-----------------------------------------------------------------------------------------
3121	c
3122	ENDIF ! fin de test sur ok_orolf
3123	C HINES GWD PARAMETRIZATION
3124
3125	IF (ok_hines) then
3126
3127	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3128	i rlat,t_seri,u_seri,v_seri,
3129	o zustrhi,zvstrhi,
3130	o d_t_hin, d_u_hin, d_v_hin)
3131	c
3132	c ajout des tendances
3133	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3134
3135	ENDIF
3136	c
3137
3138	c
3139	cIM cf. FLott BEG
3140	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3141
3142	if (mydebug) then
3143	call writefield_phy('u_seri',u_seri,llm)
3144	call writefield_phy('v_seri',v_seri,llm)
3145	call writefield_phy('t_seri',t_seri,llm)
3146	call writefield_phy('q_seri',q_seri,llm)
3147	endif
3148
3149	DO i = 1, klon
3150	zustrph(i)=0.
3151	zvstrph(i)=0.
3152	ENDDO
3153	DO k = 1, klev
3154	DO i = 1, klon
3155	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3156	c (paprs(i,k)-paprs(i,k+1))/rg
3157	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3158	c (paprs(i,k)-paprs(i,k+1))/rg
3159	ENDDO
3160	ENDDO
3161	c
3162	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3163	c
3164	IF (is_sequential) THEN
3165
3166	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3167	C ra,rg,romega,
3168	C rlat,rlon,pphis,
3169	C zustrdr,zustrli,zustrph,
3170	C zvstrdr,zvstrli,zvstrph,
3171	C paprs,u,v,
3172	C aam, torsfc)
3173	ENDIF
3174	cIM cf. FLott END
3175	cIM
3176	IF (ip_ebil_phy.ge.2) THEN
3177	ztit='after orography'
3178	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3179	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3180	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3181	END IF
3182	c
3183	c
3184	!====================================================================
3185	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3186	!====================================================================
3187	! Abderrahmane 24.08.09
3188
3189	IF (ok_cosp) THEN
3190	! adeclarer
3191	#ifdef CPP_COSP
3192	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3193
3194	print*,'freq_cosp',freq_cosp
3195	mr_ozone=wo * dobson_u * 1e3 / zmasse
3196	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3197	! s ref_liq,ref_ice
3198	call phys_cosp(itap,dtime,freq_cosp,
3199	$ ecrit_mth,ecrit_day,ecrit_hf,overlap,
3200	$ klon,klev,rlon,rlat,presnivs,
3201	$ ref_liq,ref_ice,
3202	$ pctsrf(:,is_ter)+pctsrf(:,is_lic),
3203	$ zu10m,zv10m,
3204	$ zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
3205	$ qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
3206	$ prfl(:,1:klev),psfl(:,1:klev),
3207	$ pmflxr(:,1:klev),pmflxs(:,1:klev),
3208	$ mr_ozone,cldtau, cldemi)
3209	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3210	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3211	! M clMISR,
3212	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3213	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3214
3215	ENDIF
3216
3217	#endif
3218	ENDIF !ok_cosp
3219	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3220	cAA
3221	cAA Installation de l'interface online-offline pour traceurs
3222	cAA
3223	c====================================================================
3224	c Calcul des tendances traceurs
3225	c====================================================================
3226	C
3227
3228	call phytrac (
3229	I itap, days_elapsed+1, jH_cur, debut,
3230	I lafin, dtime, u, v, t,
3231	I paprs, pplay, pmfu, pmfd,
3232	I pen_u, pde_u, pen_d, pde_d,
3233	I cdragh, coefh, fm_therm, entr_therm,
3234	I u1, v1, ftsol, pctsrf,
3235	I rlat, frac_impa, frac_nucl,rlon,
3236	I presnivs, pphis, pphi, albsol1,
3237	I qx(:,:,ivap),rhcl, cldfra, rneb,
3238	I diafra, cldliq, itop_con, ibas_con,
3239	I pmflxr, pmflxs, prfl, psfl,
3240	I da, phi, mp, upwd,
3241	I dnwd, aerosol_couple, flxmass_w,
3242	I tau_aero, piz_aero, cg_aero, ccm,
3243	I rfname,
3244	O tr_seri)
3245
3246	IF (offline) THEN
3247
3248	print*,'Attention on met a 0 les thermiques pour phystoke'
3249	call phystokenc (
3250	I nlon,klev,pdtphys,rlon,rlat,
3251	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3252	I fm_therm,entr_therm,
3253	I cdragh,coefh,u1,v1,ftsol,pctsrf,
3254	I frac_impa, frac_nucl,
3255	I pphis,airephy,dtime,itap)
3256
3257
3258	ENDIF
3259
3260	c
3261	c Calculer le transport de l'eau et de l'energie (diagnostique)
3262	c
3263	CALL transp (paprs,zxtsol,
3264	e t_seri, q_seri, u_seri, v_seri, zphi,
3265	s ve, vq, ue, uq)
3266	c
3267	cIM global posePB BEG
3268	IF(1.EQ.0) THEN
3269	c
3270	CALL transp_lay (paprs,zxtsol,
3271	e t_seri, q_seri, u_seri, v_seri, zphi,
3272	s ve_lay, vq_lay, ue_lay, uq_lay)
3273	c
3274	ENDIF !(1.EQ.0) THEN
3275	cIM global posePB END
3276	c Accumuler les variables a stocker dans les fichiers histoire:
3277	c
3278	c+jld ec_conser
3279	DO k = 1, klev
3280	DO i = 1, klon
3281	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3282	d_t_ec(i,k)=0.5/ZRCPD
3283	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3284	ENDDO
3285	ENDDO
3286
3287	DO k = 1, klev
3288	DO i = 1, klon
3289	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3290	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3291	END DO
3292	END DO
3293	c-jld ec_conser
3294	cIM
3295	IF (ip_ebil_phy.ge.1) THEN
3296	ztit='after physic'
3297	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3298	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3299	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3300	C Comme les tendances de la physique sont ajoute dans la dynamique,
3301	C on devrait avoir que la variation d'entalpie par la dynamique
3302	C est egale a la variation de la physique au pas de temps precedent.
3303	C Donc la somme de ces 2 variations devrait etre nulle.
3304
3305	call diagphy(airephy,ztit,ip_ebil_phy
3306	e , topsw, toplw, solsw, sollw, sens
3307	e , evap, rain_fall, snow_fall, ztsol
3308	e , d_h_vcol, d_qt, d_ec
3309	s , fs_bound, fq_bound )
3310	C
3311	d_h_vcol_phy=d_h_vcol
3312	C
3313	END IF
3314	C
3315	c=======================================================================
3316	c SORTIES
3317	c=======================================================================
3318
3319	cIM Interpolation sur les niveaux de pression du NMC
3320	c -------------------------------------------------
3321	c
3322	#include "calcul_STDlev.h"
3323	twriteSTD(:,:,1)=tsumSTD(:,:,2)
3324	qwriteSTD(:,:,1)=qsumSTD(:,:,2)
3325	rhwriteSTD(:,:,1)=rhsumSTD(:,:,2)
3326	phiwriteSTD(:,:,1)=phisumSTD(:,:,2)
3327	uwriteSTD(:,:,1)=usumSTD(:,:,2)
3328	vwriteSTD(:,:,1)=vsumSTD(:,:,2)
3329	wwriteSTD(:,:,1)=wsumSTD(:,:,2)
3330
3331	twriteSTD(:,:,2)=tsumSTD(:,:,1)
3332	qwriteSTD(:,:,2)=qsumSTD(:,:,1)
3333	rhwriteSTD(:,:,2)=rhsumSTD(:,:,1)
3334	phiwriteSTD(:,:,2)=phisumSTD(:,:,1)
3335	uwriteSTD(:,:,2)=usumSTD(:,:,1)
3336	vwriteSTD(:,:,2)=vsumSTD(:,:,1)
3337	wwriteSTD(:,:,2)=wsumSTD(:,:,1)
3338
3339	twriteSTD(:,:,3)=tlevSTD(:,:)
3340	qwriteSTD(:,:,3)=qlevSTD(:,:)
3341	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3342	phiwriteSTD(:,:,3)=philevSTD(:,:)
3343	uwriteSTD(:,:,3)=ulevSTD(:,:)
3344	vwriteSTD(:,:,3)=vlevSTD(:,:)
3345	wwriteSTD(:,:,3)=wlevSTD(:,:)
3346
3347	twriteSTD(:,:,4)=tlevSTD(:,:)
3348	qwriteSTD(:,:,4)=qlevSTD(:,:)
3349	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3350	phiwriteSTD(:,:,4)=philevSTD(:,:)
3351	uwriteSTD(:,:,4)=ulevSTD(:,:)
3352	vwriteSTD(:,:,4)=vlevSTD(:,:)
3353	wwriteSTD(:,:,4)=wlevSTD(:,:)
3354	c
3355	c slp sea level pressure
3356	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3357	c
3358	ccc prw = eau precipitable
3359	DO i = 1, klon
3360	prw(i) = 0.
3361	DO k = 1, klev
3362	prw(i) = prw(i) +
3363	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3364	ENDDO
3365	ENDDO
3366	c
3367	cIM initialisation + calculs divers diag AMIP2
3368	c
3369	#include "calcul_divers.h"
3370	c
3371	IF (config_inca /= 'none') THEN
3372	#ifdef INCA
3373	CALL VTe(VTphysiq)
3374	CALL VTb(VTinca)
3375
3376	CALL chemhook_end (calday,
3377	$ dtime,
3378	$ pplay,
3379	$ t_seri,
3380	$ tr_seri,
3381	$ nbtr,
3382	$ paprs,
3383	$ q_seri,
3384	$ annee_ref,
3385	$ day_ini,
3386	$ airephy,
3387	$ xjour,
3388	$ pphi,
3389	$ pphis,
3390	$ zx_rh)
3391
3392	CALL VTe(VTinca)
3393	CALL VTb(VTphysiq)
3394	#endif
3395	END IF
3396
3397	c=============================================================
3398	c
3399	c Convertir les incrementations en tendances
3400	c
3401	if (mydebug) then
3402	call writefield_phy('u_seri',u_seri,llm)
3403	call writefield_phy('v_seri',v_seri,llm)
3404	call writefield_phy('t_seri',t_seri,llm)
3405	call writefield_phy('q_seri',q_seri,llm)
3406	endif
3407
3408	DO k = 1, klev
3409	DO i = 1, klon
3410	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3411	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3412	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3413	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3414	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3415	ENDDO
3416	ENDDO
3417	c
3418	IF (nqtot.GE.3) THEN
3419	DO iq = 3, nqtot
3420	DO k = 1, klev
3421	DO i = 1, klon
3422	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3423	ENDDO
3424	ENDDO
3425	ENDDO
3426	ENDIF
3427	c
3428	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3429	cIM global posePB#include "write_bilKP_ins.h"
3430	cIM global posePB#include "write_bilKP_ave.h"
3431	c
3432	c Sauvegarder les valeurs de t et q a la fin de la physique:
3433	c
3434	DO k = 1, klev
3435	DO i = 1, klon
3436	u_ancien(i,k) = u_seri(i,k)
3437	v_ancien(i,k) = v_seri(i,k)
3438	t_ancien(i,k) = t_seri(i,k)
3439	q_ancien(i,k) = q_seri(i,k)
3440	ENDDO
3441	ENDDO
3442	c
3443	!==========================================================================
3444	! Sorties des tendances pour un point particulier
3445	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3446	! pour le debug
3447	! La valeur de igout est attribuee plus haut dans le programme
3448	!==========================================================================
3449
3450	if (prt_level.ge.1) then
3451	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3452	write(lunout,*)
3453	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3454	write(lunout,*)
3455	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
3456	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3457	s pctsrf(igout,is_sic)
3458	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3459	do k=1,klev
3460	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3461	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3462	s d_t_eva(igout,k)
3463	enddo
3464	write(lunout,*) 'cool,heat'
3465	do k=1,klev
3466	write(lunout,*) cool(igout,k),heat(igout,k)
3467	enddo
3468
3469	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3470	do k=1,klev
3471	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3472	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3473	enddo
3474
3475	write(lunout,*) 'd_ps ',d_ps(igout)
3476	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3477	do k=1,klev
3478	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3479	s d_qx(igout,k,1),d_qx(igout,k,2)
3480	enddo
3481	endif
3482
3483	!==========================================================================
3484
3485	c============================================================
3486	c Calcul de la temperature potentielle
3487	c============================================================
3488	DO k = 1, klev
3489	DO i = 1, klon
3490	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3491	ENDDO
3492	ENDDO
3493	c
3494
3495	c 22.03.04 BEG
3496	c=============================================================
3497	c Ecriture des sorties
3498	c=============================================================
3499	#ifdef CPP_IOIPSL
3500
3501	c Recupere des varibles calcule dans differents modules
3502	c pour ecriture dans histxxx.nc
3503
3504	! Get some variables from module fonte_neige_mod
3505	CALL fonte_neige_get_vars(pctsrf,
3506	. zxfqcalving, zxfqfonte, zxffonte)
3507
3508
3509	#include "phys_output_write.h"
3510
3511	#ifdef histISCCP
3512	#include "write_histISCCP.h"
3513	#endif
3514
3515	#ifdef histmthNMC
3516	#include "write_histmthNMC.h"
3517	#endif
3518
3519	#include "write_histday_seri.h"
3520
3521	#include "write_paramLMDZ_phy.h"
3522
3523	#endif
3524
3525	c 22.03.04 END
3526	c
3527	c====================================================================
3528	c Si c'est la fin, il faut conserver l'etat de redemarrage
3529	c====================================================================
3530	c
3531
3532
3533	IF (lafin) THEN
3534	itau_phy = itau_phy + itap
3535	CALL phyredem ("restartphy.nc")
3536	! open(97,form="unformatted",file="finbin")
3537	! write(97) u_seri,v_seri,t_seri,q_seri
3538	! close(97)
3539	C$OMP MASTER
3540	if (read_climoz >= 1) then
3541	if (is_mpi_root) then
3542	call nf95_close(ncid_climoz)
3543	end if
3544	deallocate(press_climoz) ! pointer
3545	end if
3546	C$OMP END MASTER
3547	ENDIF
3548
3549	! first=.false.
3550
3551	RETURN
3552	END
3553	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3554	IMPLICIT none
3555	c
3556	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3557	c la conservation de l'eau
3558	c
3559	#include "YOMCST.h"
3560	INTEGER klon,klev
3561	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3562	REAL aire(klon)
3563	REAL qtotal, zx, qcheck
3564	INTEGER i, k
3565	c
3566	zx = 0.0
3567	DO i = 1, klon
3568	zx = zx + aire(i)
3569	ENDDO
3570	qtotal = 0.0
3571	DO k = 1, klev
3572	DO i = 1, klon
3573	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3574	. *(paprs(i,k)-paprs(i,k+1))/RG
3575	ENDDO
3576	ENDDO
3577	c
3578	qcheck = qtotal/zx
3579	c
3580	RETURN
3581	END
3582	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3583	IMPLICIT none
3584	c
3585	c Tranformer une variable de la grille physique a
3586	c la grille d'ecriture
3587	c
3588	INTEGER nfield,nlon,iim,jjmp1, jjm
3589	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3590	c
3591	INTEGER i, n, ig
3592	c
3593	jjm = jjmp1 - 1
3594	DO n = 1, nfield
3595	DO i=1,iim
3596	ecrit(i,n) = fi(1,n)
3597	ecrit(i+jjm*iim,n) = fi(nlon,n)
3598	ENDDO
3599	DO ig = 1, nlon - 2
3600	ecrit(iim+ig,n) = fi(1+ig,n)
3601	ENDDO
3602	ENDDO
3603	RETURN
3604	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format