Context Navigation

physiq.F @ 1351

Last change on this file since 1351 was 1340, checked in by Laurent Fairhead, 15 years ago

Use generic rather than specific names for intrinsic procedures.
Modifications in r1334 (on Vprecip) were not propagated to physiq.F leading to a non
running program

Utilisations des noms génériques pour les procédures intrinsèques
Les modifications sur la déclaration de Vprecip dans la révision r1334 n'ont pas été
"remontées" à physiq.F d'où problème à l'exécution (en particulier sous g95)

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

File size: 117.7 KB

Line
1	! $Id: physiq.F 1340 2010-04-07 15:07:29Z musat $
2	c#define IO_DEBUG
3
4	SUBROUTINE physiq (nlon,nlev,
5	. debut,lafin,jD_cur, jH_cur,pdtphys,
6	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
7	. u,v,t,qx,
8	. flxmass_w,
9	. d_u, d_v, d_t, d_qx, d_ps
10	. , dudyn
11	. , PVteta)
12
13	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
14	$ histwrite, ju2ymds, ymds2ju, ioget_year_len
15	USE comgeomphy
16	USE phys_cal_mod
17	USE write_field_phy
18	USE dimphy
19	USE infotrac
20	USE mod_grid_phy_lmdz
21	USE mod_phys_lmdz_para
22	USE iophy
23	USE misc_mod, mydebug=>debug
24	USE vampir
25	USE pbl_surface_mod, ONLY : pbl_surface
26	USE change_srf_frac_mod
27	USE surface_data, ONLY : type_ocean, ok_veget
28	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
29	USE phys_state_var_mod ! Variables sauvegardees de la physique
30	USE phys_output_var_mod ! Variables pour les ecritures des sorties
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+1) ! 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
840	c -- convect43:
841	INTEGER ntra ! nb traceurs pour convect4.3
842	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
843	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
844	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
845	REAL dplcldt(klon), dplcldr(klon)
846	c? . condm_con(klon,klev),conda_con(klon,klev),
847	c? . mr_con(klon,klev),ep_con(klon,klev)
848	c? . ,sadiab(klon,klev),wadiab(klon,klev)
849	c --
850	c34EK
851	c
852	c Variables du changement
853	c
854	c con: convection
855	c lsc: condensation a grande echelle (Large-Scale-Condensation)
856	c ajs: ajustement sec
857	c eva: evaporation de l'eau liquide nuageuse
858	c vdf: couche limite (Vertical DiFfusion)
859	REAL rneb(klon,klev)
860
861	! tendance nulles
862	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
863
864	c
865	*********************************************************
866	* declarations
867
868	*********************************************************
869	cIM 081204 END
870	c
871	REAL pmfu(klon,klev), pmfd(klon,klev)
872	REAL pen_u(klon,klev), pen_d(klon,klev)
873	REAL pde_u(klon,klev), pde_d(klon,klev)
874	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
875	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
876	REAL prfl(klon,klev+1), psfl(klon,klev+1)
877	c
878	REAL rain_lsc(klon)
879	REAL snow_lsc(klon)
880	c
881	REAL ratqss(klon,klev),ratqsc(klon,klev)
882	real ratqsbas,ratqshaut,tau_ratqs
883	save ratqsbas,ratqshaut,tau_ratqs
884	c$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
885	real zpt_conv(klon,klev)
886
887	c Parametres lies au nouveau schema de nuages (SB, PDF)
888	real fact_cldcon
889	real facttemps
890	logical ok_newmicro
891	save ok_newmicro
892	real ref_liq(klon,klev), ref_ice(klon,klev)
893	c$OMP THREADPRIVATE(ok_newmicro)
894	save fact_cldcon,facttemps
895	c$OMP THREADPRIVATE(fact_cldcon,facttemps)
896	real facteur
897
898	integer iflag_cldcon
899	save iflag_cldcon
900	c$OMP THREADPRIVATE(iflag_cldcon)
901	logical ptconv(klon,klev)
902	cIM cf. AM 081204 BEG
903	logical ptconvth(klon,klev)
904	cIM cf. AM 081204 END
905	c
906	c Variables liees a l'ecriture de la bande histoire physique
907	c
908	c======================================================================
909	c
910	cIM cf. AM 081204 BEG
911	c declarations pour sortir sur une sous-region
912	integer imin_ins,imax_ins,jmin_ins,jmax_ins
913	save imin_ins,imax_ins,jmin_ins,jmax_ins
914	c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
915	c real lonmin_ins,lonmax_ins,latmin_ins
916	c s ,latmax_ins
917	c data lonmin_ins,lonmax_ins,latmin_ins
918	c s ,latmax_ins/
919	c valeurs initiales s -5.,20.,41.,55./
920	c s 100.,130.,-20.,20./
921	c s -180.,180.,-90.,90./
922	c======================================================================
923	cIM cf. AM 081204 END
924
925	c
926	integer itau_w ! pas de temps ecriture = itap + itau_phy
927	c
928	c
929	c Variables locales pour effectuer les appels en serie
930	c
931	REAL zx_rh(klon,klev)
932	cIM RH a 2m (la surface)
933	REAL rh2m(klon), qsat2m(klon)
934	REAL tpot(klon), tpote(klon)
935	REAL Lheat
936
937	INTEGER length
938	PARAMETER ( length = 100 )
939	REAL tabcntr0( length )
940	c
941	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
942	cIM
943	INTEGER ndex2d1(iwmax)
944	c
945	cIM AMIP2 BEG
946	REAL moyglo, mountor
947	cIM 141004 BEG
948	REAL zustrdr(klon), zvstrdr(klon)
949	REAL zustrli(klon), zvstrli(klon)
950	REAL zustrph(klon), zvstrph(klon)
951	REAL zustrhi(klon), zvstrhi(klon)
952	REAL aam, torsfc
953	cIM 141004 END
954	cIM 190504 BEG
955	INTEGER ij, imp1jmp1
956	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
957	cym A voir plus tard
958	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
959	REAL padyn(iim+1,jjmp1,klev+1)
960	REAL dudyn(iim+1,jjmp1,klev)
961	REAL rlatdyn(iim+1,jjmp1)
962	cIM 190504 END
963	LOGICAL ok_msk
964	REAL msk(klon)
965	cIM
966	REAL airetot, pi
967	cym A voir plus tard
968	cym REAL zm_wo(jjmp1, klev)
969	cIM AMIP2 END
970	c
971	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
972	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
973	c#ifdef histmthNMC
974	cym A voir plus tard !!!!
975	cym REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
976	REAL zx_tmp_fiNC(klon,nlevSTD)
977	c#endif
978	REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
979	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
980	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
981	c
982	INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
983	INTEGER nid_ctesGCM
984	SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
985	SAVE nid_ctesGCM
986	c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins, nid_nmc)
987	c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
988	c
989	cIM 280405 BEG
990	INTEGER nid_bilKPins, nid_bilKPave
991	SAVE nid_bilKPins, nid_bilKPave
992	c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
993	c
994	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
995	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
996	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
997	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
998	c
999	cIM 280405 END
1000	c
1001	INTEGER nhori, nvert, nvert1, nvert3
1002	REAL zsto, zsto1, zsto2
1003	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1004	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1005	REAL zout_isccp(napisccp)
1006	SAVE zcals, zcalh, zoutj, zout_isccp
1007	c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1008
1009	real zjulian
1010	save zjulian
1011	c$OMP THREADPRIVATE(zjulian)
1012
1013	character*20 modname
1014	character*80 abort_message
1015	logical ok_sync
1016	real date0
1017	integer idayref
1018
1019	C essai writephys
1020	integer fid_day, fid_mth, fid_ins
1021	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1022	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1023	parameter (prof2d_on = 1, prof3d_on = 2,
1024	. prof2d_av = 3, prof3d_av = 4)
1025	character*30 nom_fichier
1026	character*10 varname
1027	character*40 vartitle
1028	character*20 varunits
1029	C Variables liees au bilan d'energie et d'enthalpi
1030	REAL ztsol(klon)
1031	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1032	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1033	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1034	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1035	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1036	c$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1037	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1038	REAL d_h_vcol_phy
1039	REAL fs_bound, fq_bound
1040	SAVE d_h_vcol_phy
1041	c$OMP THREADPRIVATE(d_h_vcol_phy)
1042	REAL zero_v(klon)
1043	CHARACTER*15 ztit
1044	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
1045	SAVE ip_ebil
1046	DATA ip_ebil/0/
1047	c$OMP THREADPRIVATE(ip_ebil)
1048	INTEGER if_ebil ! level for energy conserv. dignostics
1049	SAVE if_ebil
1050	c$OMP THREADPRIVATE(if_ebil)
1051	c+jld ec_conser
1052	REAL ZRCPD
1053	c-jld ec_conser
1054	REAL t2m(klon,nbsrf) ! temperature a 2m
1055	REAL q2m(klon,nbsrf) ! humidite a 2m
1056
1057	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1058	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1059	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
1060	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1061	CHARACTER*40 tinst, tave, typeval
1062	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1063
1064	REAL re(klon, klev) ! Cloud droplet effective radius
1065	REAL fl(klon, klev) ! denominator of re
1066
1067	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1068
1069	! Aerosol optical properties
1070	CHARACTER*4, DIMENSION(naero_grp) :: rfname
1071	REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index
1072	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass concentration for all soluble aerosols[ug/m3]
1073	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi ! - " - (pre-industrial value)
1074	INTEGER :: naero ! aerosol species
1075
1076	! Parameters
1077	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1078	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1079	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1080	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
1081	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1082	! false : lecture des aerosol dans un fichier
1083	c$OMP THREADPRIVATE(aerosol_couple)
1084	INTEGER, SAVE :: flag_aerosol
1085	c$OMP THREADPRIVATE(flag_aerosol)
1086	LOGICAL, SAVE :: new_aod
1087	c$OMP THREADPRIVATE(new_aod)
1088
1089	c
1090	c Declaration des constantes et des fonctions thermodynamiques
1091	c
1092	LOGICAL,SAVE :: first=.true.
1093	c$OMP THREADPRIVATE(first)
1094
1095	integer iunit
1096
1097	integer, save:: read_climoz ! read ozone climatology
1098	C Allowed values are 0, 1 and 2
1099	C 0: do not read an ozone climatology
1100	C 1: read a single ozone climatology that will be used day and night
1101	C 2: read two ozone climatologies, the average day and night
1102	C climatology and the daylight climatology
1103
1104	integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
1105
1106	real, pointer, save:: press_climoz(:)
1107	! edges of pressure intervals for ozone climatologies, in Pa, in strictly
1108	! ascending order
1109
1110	integer, save:: co3i = 0
1111	! time index in NetCDF file of current ozone fields
1112	c$OMP THREADPRIVATE(co3i)
1113
1114	integer ro3i
1115	! required time index in NetCDF file for the ozone fields, between 1
1116	! and 360
1117
1118	#include "YOMCST.h"
1119	#include "YOETHF.h"
1120	#include "FCTTRE.h"
1121	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1122	#include "conema3.h"
1123	#include "fisrtilp.h"
1124	#include "nuage.h"
1125	#include "compbl.h"
1126	cIM 100106 END : pouvoir sortir les ctes de la physique
1127	c
1128	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1129	c Declarations pour Simulateur COSP
1130	c============================================================
1131	real :: mr_ozone(klon,klev)
1132	c======================================================================
1133	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1134	! Utilise notamment en 1D mais peut etre active egalement en 3D
1135	! en imposant la valeur de igout.
1136	c======================================================================
1137
1138	if (prt_level.ge.1) then
1139	igout=klon/2+1/klon
1140	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1141	write(lunout,*)
1142	s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1143	write(lunout,*)
1144	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1145
1146	write(lunout,*) 'paprs, play, phi, u, v, t'
1147	do k=1,klev
1148	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1149	s u(igout,k),v(igout,k),t(igout,k)
1150	enddo
1151	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1152	do k=1,klev
1153	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1154	enddo
1155	endif
1156
1157	c======================================================================
1158
1159	cym => necessaire pour iflag_con != 2
1160	pmfd(:,:) = 0.
1161	pen_u(:,:) = 0.
1162	pen_d(:,:) = 0.
1163	pde_d(:,:) = 0.
1164	pde_u(:,:) = 0.
1165	aam=0.
1166
1167	torsfc=0.
1168	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1169
1170	if (first) then
1171
1172	cCR:nvelles variables convection/poches froides
1173
1174	print*, '================================================='
1175	print*, 'Allocation des variables locales et sauvegardees'
1176	call phys_local_var_init
1177	c appel a la lecture du run.def physique
1178	call conf_phys(ok_journe, ok_mensuel,
1179	. ok_instan, ok_hf,
1180	. ok_LES,
1181	. solarlong0,seuil_inversion,
1182	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1183	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
1184	. ok_ade, ok_aie, aerosol_couple,
1185	. flag_aerosol, new_aod,
1186	. bl95_b0, bl95_b1,
1187	. iflag_thermals,nsplit_thermals,tau_thermals,
1188	. iflag_thermals_ed,iflag_thermals_optflux,
1189	c nv flags pour la convection et les poches froides
1190	. iflag_coupl,iflag_clos,iflag_wake, read_climoz)
1191	call phys_state_var_init(read_climoz)
1192	call phys_output_var_init
1193	print*, '================================================='
1194
1195	cIM beg
1196	dnwd0=0.0
1197	ftd=0.0
1198	fqd=0.0
1199	cin=0.
1200	cym Attention pbase pas initialise dans concvl !!!!
1201	pbase=0
1202	paire_ter(:)=0.
1203	cIM 180608
1204	c pmflxr=0.
1205	c pmflxs=0.
1206	itau_con=0
1207	first=.false.
1208
1209	endif ! first
1210
1211	modname = 'physiq'
1212	cIM
1213	IF (ip_ebil_phy.ge.1) THEN
1214	DO i=1,klon
1215	zero_v(i)=0.
1216	END DO
1217	END IF
1218	ok_sync=.TRUE.
1219
1220	IF (debut) THEN
1221	CALL suphel ! initialiser constantes et parametres phys.
1222	ENDIF
1223
1224	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1225
1226
1227	c======================================================================
1228	! Gestion calendrier : mise a jour du module phys_cal_mod
1229	!
1230	CALL phys_cal_update(jD_cur,jH_cur)
1231
1232	c
1233	c Si c'est le debut, il faut initialiser plusieurs choses
1234	c ********
1235	c
1236	IF (debut) THEN
1237	!rv
1238	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1239	cde la convection a partir des caracteristiques du thermique
1240	wght_th(:,:)=1.
1241	lalim_conv(:)=1
1242	cRC
1243	u10m(:,:)=0.
1244	v10m(:,:)=0.
1245	rain_con(:)=0.
1246	snow_con(:)=0.
1247	topswai(:)=0.
1248	topswad(:)=0.
1249	solswai(:)=0.
1250	solswad(:)=0.
1251
1252	lambda_th(:,:)=0.
1253	wmax_th(:)=0.
1254	tau_overturning_th(:)=0.
1255
1256	IF (config_inca /= 'none') THEN
1257	! jg : initialisation jusqu'au ces variables sont dans restart
1258	ccm(:,:,:) = 0.
1259	tau_aero(:,:,:,:) = 0.
1260	piz_aero(:,:,:,:) = 0.
1261	cg_aero(:,:,:,:) = 0.
1262	END IF
1263
1264	rnebcon0(:,:) = 0.0
1265	clwcon0(:,:) = 0.0
1266	rnebcon(:,:) = 0.0
1267	clwcon(:,:) = 0.0
1268
1269	cIM
1270	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1271	c
1272	print*,'iflag_coupl,iflag_clos,iflag_wake',
1273	. iflag_coupl,iflag_clos,iflag_wake
1274	print*,'CYCLE_DIURNE', cycle_diurne
1275	c
1276	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1277	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1278	CALL abort_gcm (modname,abort_message,1)
1279	ENDIF
1280	c
1281	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1282	abort_message = 'ISCCP-like outputs may be available for KE
1283	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1284	CALL abort_gcm (modname,abort_message,1)
1285	ENDIF
1286	c
1287	c Initialiser les compteurs:
1288	c
1289	itap = 0
1290	itaprad = 0
1291
1292	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1293	!! Un petit travail à faire ici.
1294	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1295
1296	if (iflag_pbl>1) then
1297	PRINT*, "Using method MELLOR&YAMADA"
1298	endif
1299
1300	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1301	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1302	! dyn3d
1303	! Attention : la version precedente n'etait pas tres propre.
1304	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1305	! pour obtenir le meme resultat.
1306	dtime=pdtphys
1307	radpas = NINT( 86400./dtime/nbapp_rad)
1308	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1309
1310	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1311	cIM begin
1312	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1313	$,ratqs(1,1)
1314	cIM end
1315
1316
1317
1318	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1319	c
1320	C on remet le calendrier a zero
1321	c
1322	IF (raz_date .eq. 1) THEN
1323	itau_phy = 0
1324	ENDIF
1325
1326	cIM cf. AM 081204 BEG
1327	PRINT*,'cycle_diurne3 =',cycle_diurne
1328	cIM cf. AM 081204 END
1329	c
1330	CALL printflag( tabcntr0,radpas,ok_journe,
1331	, ok_instan, ok_region )
1332	c
1333	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1334	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1335	. pdtphys
1336	abort_message='Pas physique n est pas correct '
1337	! call abort_gcm(modname,abort_message,1)
1338	dtime=pdtphys
1339	ENDIF
1340	IF (nlon .NE. klon) THEN
1341	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1342	. klon
1343	abort_message='nlon et klon ne sont pas coherents'
1344	call abort_gcm(modname,abort_message,1)
1345	ENDIF
1346	IF (nlev .NE. klev) THEN
1347	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1348	. klev
1349	abort_message='nlev et klev ne sont pas coherents'
1350	call abort_gcm(modname,abort_message,1)
1351	ENDIF
1352	c
1353	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1354	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1355	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1356	abort_message='Nbre d appels au rayonnement insuffisant'
1357	call abort_gcm(modname,abort_message,1)
1358	ENDIF
1359	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1360	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1361	. ok_cvl
1362	c
1363	cKE43
1364	c Initialisation pour la convection de K.E. (sb):
1365	IF (iflag_con.GE.3) THEN
1366
1367	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1368	WRITE(lunout,*)
1369	. "On va utiliser le melange convectif des traceurs qui"
1370	WRITE(lunout,*)"est calcule dans convect4.3"
1371	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1372
1373	DO i = 1, klon
1374	ema_cbmf(i) = 0.
1375	ema_pcb(i) = 0.
1376	ema_pct(i) = 0.
1377	ema_workcbmf(i) = 0.
1378	ENDDO
1379	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1380	DO i = 1, klon
1381	ibas_con(i) = 1
1382	itop_con(i) = 1
1383	ENDDO
1384	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1385	c===============================================================================
1386	cCR:04.12.07: initialisations poches froides
1387	c Controle de ALE et ALP pour la fermeture convective (jyg)
1388	if (iflag_wake.eq.1) then
1389	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1390	s ,alp_bl_prescr, ale_bl_prescr)
1391	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1392	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1393	endif
1394
1395	do i = 1,klon
1396	Ale_bl(i)=0.
1397	Alp_bl(i)=0.
1398	enddo
1399
1400	c================================================================================
1401
1402	ENDIF
1403
1404	DO i=1,klon
1405	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1406	ENDDO
1407
1408	c34EK
1409	IF (ok_orodr) THEN
1410
1411	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1412	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1413	! justement quand ok_orodr = false.
1414	! ce rugoro est utilise par la couche limite et fait double emploi
1415	! avec les paramétrisations spécifiques de Francois Lott.
1416	! DO i=1,klon
1417	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1418	! ENDDO
1419	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1420	IF (ok_strato) THEN
1421	CALL SUGWD_strato(klon,klev,paprs,pplay)
1422	ELSE
1423	CALL SUGWD(klon,klev,paprs,pplay)
1424	ENDIF
1425
1426	DO i=1,klon
1427	zuthe(i)=0.
1428	zvthe(i)=0.
1429	if(zstd(i).gt.10.)then
1430	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1431	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1432	endif
1433	ENDDO
1434	ENDIF
1435	c
1436	c
1437	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1438	WRITE(lunout,*)'La frequence de lecture surface est de ',
1439	. lmt_pas
1440	c
1441	cIM 030306 END
1442
1443	capemaxcels = 't_max(X)'
1444	t2mincels = 't_min(X)'
1445	t2maxcels = 't_max(X)'
1446	tinst = 'inst(X)'
1447	tave = 'ave(X)'
1448	cIM cf. AM 081204 BEG
1449	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1450	cIM cf. AM 081204 END
1451	c
1452	c=============================================================
1453	c Initialisation des sorties
1454	c=============================================================
1455
1456	#ifdef CPP_IOIPSL
1457
1458	c$OMP MASTER
1459	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1460	& ctetaSTD,dtime,ok_veget,
1461	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1462	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
1463	& read_climoz, new_aod, aerosol_couple)
1464	c$OMP END MASTER
1465	c$OMP BARRIER
1466
1467	#ifdef histISCCP
1468	#include "ini_histISCCP.h"
1469	#endif
1470
1471	#ifdef histmthNMC
1472	#include "ini_histmthNMC.h"
1473	#endif
1474
1475	#include "ini_histday_seri.h"
1476
1477	#include "ini_paramLMDZ_phy.h"
1478
1479	#endif
1480
1481	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1482	ecrit_hf2mth = ecrit_mth/ecrit_hf
1483
1484	ecrit_hf = ecrit_hf * un_jour
1485	!IM
1486	IF(ecrit_day.LE.1.) THEN
1487	ecrit_day = ecrit_day * un_jour !en secondes
1488	ENDIF
1489	!IM
1490	ecrit_mth = ecrit_mth * un_jour
1491	ecrit_ins = ecrit_ins * un_jour
1492	ecrit_reg = ecrit_reg * un_jour
1493	ecrit_tra = ecrit_tra * un_jour
1494	ecrit_ISCCP = ecrit_ISCCP * un_jour
1495	ecrit_LES = ecrit_LES * un_jour
1496	c
1497	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1498	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1499	. ecrit_hf2mth
1500	cIM 030306 END
1501
1502
1503	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1504	date0 = jD_ref
1505	WRITE(,) 'physiq date0 : ',date0
1506	c
1507	c
1508	c
1509	c Prescrire l'ozone dans l'atmosphere
1510	c
1511	c
1512	cc DO i = 1, klon
1513	cc DO k = 1, klev
1514	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1515	cc ENDDO
1516	cc ENDDO
1517	c
1518	IF (config_inca /= 'none') THEN
1519	#ifdef INCA
1520	CALL VTe(VTphysiq)
1521	CALL VTb(VTinca)
1522	! iii = MOD(NINT(xjour),360)
1523	! calday = FLOAT(iii) + jH_cur
1524	calday = FLOAT(days_elapsed) + jH_cur
1525	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1526
1527	CALL chemini(
1528	$ rg,
1529	$ ra,
1530	$ airephy,
1531	$ rlat,
1532	$ rlon,
1533	$ presnivs,
1534	$ calday,
1535	$ klon,
1536	$ nqtot,
1537	$ pdtphys,
1538	$ annee_ref,
1539	$ day_ref,
1540	$ itau_phy)
1541
1542	CALL VTe(VTinca)
1543	CALL VTb(VTphysiq)
1544	#endif
1545	END IF
1546	c
1547	c
1548	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1549	! Nouvelle initialisation pour le rayonnement RRTM
1550	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1551
1552	call iniradia(klon,klev,paprs(1,1:klev+1))
1553
1554	C$omp single
1555	if (read_climoz >= 1) then
1556	call open_climoz(ncid_climoz, press_climoz)
1557	END IF
1558	C$omp end single
1559	ENDIF
1560	!
1561	! ************** Fin de IF ( debut ) *************
1562	!
1563	!
1564	! Incrementer le compteur de la physique
1565	!
1566	itap = itap + 1
1567	!
1568	! Update fraction of the sub-surfaces (pctsrf) and
1569	! initialize, where a new fraction has appeared, all variables depending
1570	! on the surface fraction.
1571	!
1572	CALL change_srf_frac(itap, dtime, days_elapsed+1,
1573	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1574
1575	! Tendances bidons pour les processus qui n'affectent pas certaines
1576	! variables.
1577	du0(:,:)=0.
1578	dv0(:,:)=0.
1579	dq0(:,:)=0.
1580	dql0(:,:)=0.
1581	c
1582	c Mettre a zero des variables de sortie (pour securite)
1583	c
1584	DO i = 1, klon
1585	d_ps(i) = 0.0
1586	ENDDO
1587	DO k = 1, klev
1588	DO i = 1, klon
1589	d_t(i,k) = 0.0
1590	d_u(i,k) = 0.0
1591	d_v(i,k) = 0.0
1592	ENDDO
1593	ENDDO
1594	DO iq = 1, nqtot
1595	DO k = 1, klev
1596	DO i = 1, klon
1597	d_qx(i,k,iq) = 0.0
1598	ENDDO
1599	ENDDO
1600	ENDDO
1601	da(:,:)=0.
1602	mp(:,:)=0.
1603	phi(:,:,:)=0.
1604	c
1605	c Ne pas affecter les valeurs entrees de u, v, h, et q
1606	c
1607	DO k = 1, klev
1608	DO i = 1, klon
1609	t_seri(i,k) = t(i,k)
1610	u_seri(i,k) = u(i,k)
1611	v_seri(i,k) = v(i,k)
1612	q_seri(i,k) = qx(i,k,ivap)
1613	ql_seri(i,k) = qx(i,k,iliq)
1614	qs_seri(i,k) = 0.
1615	ENDDO
1616	ENDDO
1617	IF (nqtot.GE.3) THEN
1618	DO iq = 3, nqtot
1619	DO k = 1, klev
1620	DO i = 1, klon
1621	tr_seri(i,k,iq-2) = qx(i,k,iq)
1622	ENDDO
1623	ENDDO
1624	ENDDO
1625	ELSE
1626	DO k = 1, klev
1627	DO i = 1, klon
1628	tr_seri(i,k,1) = 0.0
1629	ENDDO
1630	ENDDO
1631	ENDIF
1632	C
1633	DO i = 1, klon
1634	ztsol(i) = 0.
1635	ENDDO
1636	DO nsrf = 1, nbsrf
1637	DO i = 1, klon
1638	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1639	ENDDO
1640	ENDDO
1641	cIM
1642	IF (ip_ebil_phy.ge.1) THEN
1643	ztit='after dynamic'
1644	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1645	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1646	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1647	C Comme les tendances de la physique sont ajoute dans la dynamique,
1648	C on devrait avoir que la variation d'entalpie par la dynamique
1649	C est egale a la variation de la physique au pas de temps precedent.
1650	C Donc la somme de ces 2 variations devrait etre nulle.
1651	call diagphy(airephy,ztit,ip_ebil_phy
1652	e , zero_v, zero_v, zero_v, zero_v, zero_v
1653	e , zero_v, zero_v, zero_v, ztsol
1654	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1655	s , fs_bound, fq_bound )
1656	END IF
1657
1658	c Diagnostiquer la tendance dynamique
1659	c
1660	IF (ancien_ok) THEN
1661	DO k = 1, klev
1662	DO i = 1, klon
1663	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1664	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1665	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1666	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1667	ENDDO
1668	ENDDO
1669	ELSE
1670	DO k = 1, klev
1671	DO i = 1, klon
1672	d_u_dyn(i,k) = 0.0
1673	d_v_dyn(i,k) = 0.0
1674	d_t_dyn(i,k) = 0.0
1675	d_q_dyn(i,k) = 0.0
1676	ENDDO
1677	ENDDO
1678	ancien_ok = .TRUE.
1679	ENDIF
1680	c
1681	c Ajouter le geopotentiel du sol:
1682	c
1683	DO k = 1, klev
1684	DO i = 1, klon
1685	zphi(i,k) = pphi(i,k) + pphis(i)
1686	ENDDO
1687	ENDDO
1688	c
1689	c Verifier les temperatures
1690	c
1691	cIM BEG
1692	IF (check) THEN
1693	amn=MIN(ftsol(1,is_ter),1000.)
1694	amx=MAX(ftsol(1,is_ter),-1000.)
1695	DO i=2, klon
1696	amn=MIN(ftsol(i,is_ter),amn)
1697	amx=MAX(ftsol(i,is_ter),amx)
1698	ENDDO
1699	c
1700	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1701	ENDIF !(check) THEN
1702	cIM END
1703	c
1704	CALL hgardfou(t_seri,ftsol,'debutphy')
1705	c
1706	cIM BEG
1707	IF (check) THEN
1708	amn=MIN(ftsol(1,is_ter),1000.)
1709	amx=MAX(ftsol(1,is_ter),-1000.)
1710	DO i=2, klon
1711	amn=MIN(ftsol(i,is_ter),amn)
1712	amx=MAX(ftsol(i,is_ter),amx)
1713	ENDDO
1714	c
1715	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1716	ENDIF !(check) THEN
1717	cIM END
1718	c
1719	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1720	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1721	c
1722	if (read_climoz >= 1) then
1723	C Ozone from a file
1724	! Update required ozone index:
1725	ro3i = int((days_elapsed + jh_cur - jh_1jan)
1726	$ / ioget_year_len(year_cur) * 360.) + 1
1727	if (ro3i == 361) ro3i = 360
1728	C (This should never occur, except perhaps because of roundup
1729	C error. See documentation.)
1730	if (ro3i /= co3i) then
1731	C Update ozone field:
1732	if (read_climoz == 1) then
1733	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1734	$ press_in_edg=press_climoz, paprs=paprs, v3=wo)
1735	else
1736	C read_climoz == 2
1737	call regr_pr_av(ncid_climoz,
1738	$ (/"tro3 ", "tro3_daylight"/),
1739	$ julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1740	$ v3=wo)
1741	end if
1742	! Convert from mole fraction of ozone to column density of ozone in a
1743	! cell, in kDU:
1744	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1745	$ * rmo3 / rmd * zmasse / dobson_u / 1e3
1746	C (By regridding ozone values for LMDZ only once every 360th of
1747	C year, we have already neglected the variation of pressure in one
1748	C 360th of year. So do not recompute "wo" at each time step even if
1749	C "zmasse" changes a little.)
1750	co3i = ro3i
1751	end if
1752	elseif (MOD(itap-1,lmt_pas) == 0) THEN
1753	C Once per day, update ozone from Royer:
1754	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1755	ENDIF
1756	c
1757	c Re-evaporer l'eau liquide nuageuse
1758	c
1759	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1760	DO i = 1, klon
1761	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1762	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1763	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1764	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1765	zb = MAX(0.0,ql_seri(i,k))
1766	za = - MAX(0.0,ql_seri(i,k))
1767	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1768	t_seri(i,k) = t_seri(i,k) + za
1769	q_seri(i,k) = q_seri(i,k) + zb
1770	ql_seri(i,k) = 0.0
1771	d_t_eva(i,k) = za
1772	d_q_eva(i,k) = zb
1773	ENDDO
1774	ENDDO
1775	cIM
1776	IF (ip_ebil_phy.ge.2) THEN
1777	ztit='after reevap'
1778	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1779	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1780	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1781	call diagphy(airephy,ztit,ip_ebil_phy
1782	e , zero_v, zero_v, zero_v, zero_v, zero_v
1783	e , zero_v, zero_v, zero_v, ztsol
1784	e , d_h_vcol, d_qt, d_ec
1785	s , fs_bound, fq_bound )
1786	C
1787	END IF
1788
1789	c
1790	c=========================================================================
1791	! Calculs de l'orbite.
1792	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1793	! doit donc etre placé avant radlwsw et pbl_surface
1794
1795	! calcul selon la routine utilisee pour les planetes
1796	if (new_orbit) then
1797	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1798	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1799	! day_since_equinox = (jD_cur) - jD_eq
1800	call solarlong(day_since_equinox, zlongi, dist)
1801	else
1802	! calcul selon la routine utilisee pour l'AR4
1803	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1804	! solarlong0
1805	if (solarlong0<-999.) then
1806	CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
1807	else
1808	zlongi=solarlong0 ! longitude solaire vraie
1809	dist=1. ! distance au soleil / moyenne
1810	endif
1811	endif
1812	if(prt_level.ge.1) &
1813	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1814
1815	! Avec ou sans cycle diurne
1816	IF (cycle_diurne) THEN
1817	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1818	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1819	ELSE
1820	CALL angle(zlongi, rlat, fract, rmu0)
1821	ENDIF
1822
1823	if (mydebug) then
1824	call writefield_phy('u_seri',u_seri,llm)
1825	call writefield_phy('v_seri',v_seri,llm)
1826	call writefield_phy('t_seri',t_seri,llm)
1827	call writefield_phy('q_seri',q_seri,llm)
1828	endif
1829
1830	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1831	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1832	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1833	c turbulent du couche limit.
1834	c
1835	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1836	c ecriture des fihiers hist_XXXX.nc, ces sont :
1837	c qsol, zq2m, s_pblh, s_lcl,
1838	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1839	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1840	c zxrugs, zu10m, zv10m, fder,
1841	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1842	c frugs, agesno, fsollw, fsolsw,
1843	c d_ts, fevap, fluxlat, t2m,
1844	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1845	c
1846	c Certains ne sont pas utiliser du tout :
1847	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1848	c
1849
1850	CALL pbl_surface(
1851	e dtime, date0, itap, days_elapsed+1,
1852	e debut, lafin,
1853	e rlon, rlat, rugoro, rmu0,
1854	e rain_fall, snow_fall, solsw, sollw,
1855	e t_seri, q_seri, u_seri, v_seri,
1856	e pplay, paprs, pctsrf,
1857	+ ftsol, falb1, falb2, u10m, v10m,
1858	s sollwdown, cdragh, cdragm, u1, v1,
1859	s albsol1, albsol2, sens, evap,
1860	s zxtsol, zxfluxlat, zt2m, qsat2m,
1861	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1862	s coefh, slab_wfbils,
1863	d qsol, zq2m, s_pblh, s_lcl,
1864	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1865	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1866	d zxrugs, zu10m, zv10m, fder,
1867	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1868	d frugs, agesno, fsollw, fsolsw,
1869	d d_ts, fevap, fluxlat, t2m,
1870	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1871	- dsens, devap, zxsnow,
1872	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1873
1874
1875	!-----------------------------------------------------------------------------------------
1876	! ajout des tendances de la diffusion turbulente
1877	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1878	!-----------------------------------------------------------------------------------------
1879
1880	if (mydebug) then
1881	call writefield_phy('u_seri',u_seri,llm)
1882	call writefield_phy('v_seri',v_seri,llm)
1883	call writefield_phy('t_seri',t_seri,llm)
1884	call writefield_phy('q_seri',q_seri,llm)
1885	endif
1886
1887
1888	IF (ip_ebil_phy.ge.2) THEN
1889	ztit='after surface_main'
1890	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1891	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1892	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1893	call diagphy(airephy,ztit,ip_ebil_phy
1894	e , zero_v, zero_v, zero_v, zero_v, sens
1895	e , evap , zero_v, zero_v, ztsol
1896	e , d_h_vcol, d_qt, d_ec
1897	s , fs_bound, fq_bound )
1898	END IF
1899
1900	c =================================================================== c
1901	c Calcul de Qsat
1902
1903	DO k = 1, klev
1904	DO i = 1, klon
1905	zx_t = t_seri(i,k)
1906	IF (thermcep) THEN
1907	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1908	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1909	zx_qs = MIN(0.5,zx_qs)
1910	zcor = 1./(1.-retv*zx_qs)
1911	zx_qs = zx_qs*zcor
1912	ELSE
1913	IF (zx_t.LT.t_coup) THEN
1914	zx_qs = qsats(zx_t)/pplay(i,k)
1915	ELSE
1916	zx_qs = qsatl(zx_t)/pplay(i,k)
1917	ENDIF
1918	ENDIF
1919	zqsat(i,k)=zx_qs
1920	ENDDO
1921	ENDDO
1922
1923	if (prt_level.ge.1) then
1924	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1925	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1926	endif
1927	c
1928	c Appeler la convection (au choix)
1929	c
1930	DO k = 1, klev
1931	DO i = 1, klon
1932	conv_q(i,k) = d_q_dyn(i,k)
1933	. + d_q_vdf(i,k)/dtime
1934	conv_t(i,k) = d_t_dyn(i,k)
1935	. + d_t_vdf(i,k)/dtime
1936	ENDDO
1937	ENDDO
1938	IF (check) THEN
1939	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1940	WRITE(lunout,*) "avantcon=", za
1941	ENDIF
1942	zx_ajustq = .FALSE.
1943	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1944	IF (zx_ajustq) THEN
1945	DO i = 1, klon
1946	z_avant(i) = 0.0
1947	ENDDO
1948	DO k = 1, klev
1949	DO i = 1, klon
1950	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1951	. *(paprs(i,k)-paprs(i,k+1))/RG
1952	ENDDO
1953	ENDDO
1954	ENDIF
1955
1956	c Calcule de vitesse verticale a partir de flux de masse verticale
1957	DO k = 1, klev
1958	DO i = 1, klon
1959	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1960	END DO
1961	END DO
1962	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
1963	$ omega(igout, :)
1964
1965	IF (iflag_con.EQ.1) THEN
1966	stop'reactiver le call conlmd dans physiq.F'
1967	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1968	c . d_t_con, d_q_con,
1969	c . rain_con, snow_con, ibas_con, itop_con)
1970	ELSE IF (iflag_con.EQ.2) THEN
1971	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1972	e conv_t, conv_q, -evap, omega,
1973	s d_t_con, d_q_con, rain_con, snow_con,
1974	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1975	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1976	d_u_con = 0.
1977	d_v_con = 0.
1978
1979	WHERE (rain_con < 0.) rain_con = 0.
1980	WHERE (snow_con < 0.) snow_con = 0.
1981	DO i = 1, klon
1982	ibas_con(i) = klev+1 - kcbot(i)
1983	itop_con(i) = klev+1 - kctop(i)
1984	ENDDO
1985	ELSE IF (iflag_con.GE.3) THEN
1986	c nb of tracers for the KE convection:
1987	c MAF la partie traceurs est faite dans phytrac
1988	c on met ntra=1 pour limiter les appels mais on peut
1989	c supprimer les calculs / ftra.
1990	ntra = 1
1991
1992	c=====================================================================================
1993	cajout pour la parametrisation des poches froides:
1994	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1995	do k=1,klev
1996	do i=1,klon
1997	if (iflag_wake.eq.1) then
1998	t_wake(i,k) = t_seri(i,k)
1999	. +(1-wake_s(i))*wake_deltat(i,k)
2000	q_wake(i,k) = q_seri(i,k)
2001	. +(1-wake_s(i))*wake_deltaq(i,k)
2002	t_undi(i,k) = t_seri(i,k)
2003	. -wake_s(i)*wake_deltat(i,k)
2004	q_undi(i,k) = q_seri(i,k)
2005	. -wake_s(i)*wake_deltaq(i,k)
2006	else
2007	t_wake(i,k) = t_seri(i,k)
2008	q_wake(i,k) = q_seri(i,k)
2009	t_undi(i,k) = t_seri(i,k)
2010	q_undi(i,k) = q_seri(i,k)
2011	endif
2012	enddo
2013	enddo
2014
2015	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2016	cc-- pour le soulevement des particules dans le modele convectif
2017	c
2018	do i = 1,klon
2019	ALE(i) = 0.
2020	ALP(i) = 0.
2021	enddo
2022	c
2023	ccalcul de ale_wake et alp_wake
2024	do i = 1,klon
2025	if (iflag_wake.eq.1) then
2026	ale_wake(i) = 0.5wake_cstar(i)*2
2027	alp_wake(i) = wake_fip(i)
2028	else
2029	ale_wake(i) = 0.
2030	alp_wake(i) = 0.
2031	endif
2032	enddo
2033	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2034	cdans le thermique sinon
2035	if (iflag_coupl.eq.0) then
2036	if (debut) print*,'ALE et ALP imposes'
2037	do i = 1,klon
2038	con ne couple que ale
2039	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2040	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2041	con ne couple que alp
2042	c ALP(i) = alp_wake(i) + Alp_bl(i)
2043	ALP(i) = alp_wake(i) + alp_bl_prescr
2044	enddo
2045	else
2046	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2047	do i = 1,klon
2048	ALE(i) = max(ale_wake(i),Ale_bl(i))
2049	ALP(i) = alp_wake(i) + Alp_bl(i)
2050	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2051	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2052	enddo
2053	endif
2054	do i=1,klon
2055	if (alp(i)>alp_max) then
2056	IF(prt_level>9)WRITE(lunout,*) &
2057	& 'WARNING SUPER ALP (seuil=',alp_max,
2058	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2059	alp(i)=alp_max
2060	endif
2061	if (ale(i)>ale_max) then
2062	IF(prt_level>9)WRITE(lunout,*) &
2063	& 'WARNING SUPER ALE (seuil=',ale_max,
2064	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2065	ale(i)=ale_max
2066	endif
2067	enddo
2068
2069	cfin calcul ale et alp
2070	c=================================================================================================
2071
2072
2073	c sb, oct02:
2074	c Schema de convection modularise et vectorise:
2075	c (driver commun aux versions 3 et 4)
2076	c
2077	IF (ok_cvl) THEN ! new driver for convectL
2078
2079	CALL concvl (iflag_con,iflag_clos,
2080	. dtime,paprs,pplay,t_undi,q_undi,
2081	. t_wake,q_wake,wake_s,
2082	. u_seri,v_seri,tr_seri,nbtr,
2083	. ALE,ALP,
2084	. ema_work1,ema_work2,
2085	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2086	. rain_con, snow_con, ibas_con, itop_con, sigd,
2087	. upwd,dnwd,dnwd0,
2088	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2089	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2090	. pmflxr,pmflxs,da,phi,mp,
2091	. ftd,fqd,lalim_conv,wght_th)
2092
2093	cIM begin
2094	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2095	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2096	cIM end
2097	cIM cf. FH
2098	clwcon0=qcondc
2099	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2100
2101	do i = 1, klon
2102	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2103	enddo
2104
2105	ELSE ! ok_cvl
2106	c MAF conema3 ne contient pas les traceurs
2107	CALL conema3 (dtime,
2108	. paprs,pplay,t_seri,q_seri,
2109	. u_seri,v_seri,tr_seri,ntra,
2110	. ema_work1,ema_work2,
2111	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2112	. rain_con, snow_con, ibas_con, itop_con,
2113	. upwd,dnwd,dnwd0,bas,top,
2114	. Ma,cape,tvp,rflag,
2115	. pbase
2116	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2117	. ,clwcon0)
2118
2119	ENDIF ! ok_cvl
2120
2121	c
2122	c Correction precip
2123	rain_con = rain_con * cvl_corr
2124	snow_con = snow_con * cvl_corr
2125	c
2126
2127	IF (.NOT. ok_gust) THEN
2128	do i = 1, klon
2129	wd(i)=0.0
2130	enddo
2131	ENDIF
2132
2133	c =================================================================== c
2134	c Calcul des proprietes des nuages convectifs
2135	c
2136
2137	c calcul des proprietes des nuages convectifs
2138	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2139	call clouds_gno
2140	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2141
2142	c =================================================================== c
2143
2144	DO i = 1, klon
2145	ema_pcb(i) = paprs(i,ibas_con(i))
2146	ENDDO
2147	DO i = 1, klon
2148	! L'idicage de itop_con peut cacher un pb potentiel
2149	! FH sous la dictee de JYG, CR
2150	ema_pct(i) = paprs(i,itop_con(i)+1)
2151
2152	if (itop_con(i).gt.klev-3) then
2153	if(prt_level >= 9) then
2154	write(lunout,*)'La convection monte trop haut '
2155	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2156	endif
2157	endif
2158	ENDDO
2159	DO i = 1, klon
2160	ema_cbmf(i) = ema_workcbmf(i)
2161	ENDDO
2162	ELSE IF (iflag_con.eq.0) THEN
2163	write(lunout,*) 'On n appelle pas la convection'
2164	clwcon0=0.
2165	rnebcon0=0.
2166	d_t_con=0.
2167	d_q_con=0.
2168	d_u_con=0.
2169	d_v_con=0.
2170	rain_con=0.
2171	snow_con=0.
2172	bas=1
2173	top=1
2174	ELSE
2175	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2176	CALL abort
2177	ENDIF
2178
2179	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2180	c . d_u_con, d_v_con)
2181
2182	!-----------------------------------------------------------------------------------------
2183	! ajout des tendances de la diffusion turbulente
2184	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2185	!-----------------------------------------------------------------------------------------
2186
2187	if (mydebug) then
2188	call writefield_phy('u_seri',u_seri,llm)
2189	call writefield_phy('v_seri',v_seri,llm)
2190	call writefield_phy('t_seri',t_seri,llm)
2191	call writefield_phy('q_seri',q_seri,llm)
2192	endif
2193
2194	cIM
2195	IF (ip_ebil_phy.ge.2) THEN
2196	ztit='after convect'
2197	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2198	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2199	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2200	call diagphy(airephy,ztit,ip_ebil_phy
2201	e , zero_v, zero_v, zero_v, zero_v, zero_v
2202	e , zero_v, rain_con, snow_con, ztsol
2203	e , d_h_vcol, d_qt, d_ec
2204	s , fs_bound, fq_bound )
2205	END IF
2206	C
2207	IF (check) THEN
2208	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2209	WRITE(lunout,*)"aprescon=", za
2210	zx_t = 0.0
2211	za = 0.0
2212	DO i = 1, klon
2213	za = za + airephy(i)/FLOAT(klon)
2214	zx_t = zx_t + (rain_con(i)+
2215	. snow_con(i))*airephy(i)/FLOAT(klon)
2216	ENDDO
2217	zx_t = zx_t/za*dtime
2218	WRITE(lunout,*)"Precip=", zx_t
2219	ENDIF
2220	IF (zx_ajustq) THEN
2221	DO i = 1, klon
2222	z_apres(i) = 0.0
2223	ENDDO
2224	DO k = 1, klev
2225	DO i = 1, klon
2226	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2227	. *(paprs(i,k)-paprs(i,k+1))/RG
2228	ENDDO
2229	ENDDO
2230	DO i = 1, klon
2231	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2232	. /z_apres(i)
2233	ENDDO
2234	DO k = 1, klev
2235	DO i = 1, klon
2236	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2237	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2238	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2239	ENDIF
2240	ENDDO
2241	ENDDO
2242	ENDIF
2243	zx_ajustq=.FALSE.
2244
2245	c
2246	c=============================================================================
2247	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2248	cpour la couche limite diffuse pour l instant
2249	c
2250	if (iflag_wake.eq.1) then
2251	DO k=1,klev
2252	DO i=1,klon
2253	dt_dwn(i,k) = ftd(i,k)
2254	wdt_PBL(i,k) = 0.
2255	dq_dwn(i,k) = fqd(i,k)
2256	wdq_PBL(i,k) = 0.
2257	M_dwn(i,k) = dnwd0(i,k)
2258	M_up(i,k) = upwd(i,k)
2259	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2260	udt_PBL(i,k) = 0.
2261	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2262	udq_PBL(i,k) = 0.
2263	ENDDO
2264	ENDDO
2265	c
2266	ccalcul caracteristiques de la poche froide
2267	call calWAKE (paprs,pplay,dtime
2268	: ,t_seri,q_seri,omega
2269	: ,dt_dwn,dq_dwn,M_dwn,M_up
2270	: ,dt_a,dq_a,sigd
2271	: ,wdt_PBL,wdq_PBL
2272	: ,udt_PBL,udq_PBL
2273	o ,wake_deltat,wake_deltaq,wake_dth
2274	o ,wake_h,wake_s,wake_dens
2275	o ,wake_pe,wake_fip,wake_gfl
2276	o ,dt_wake,dq_wake
2277	o ,wake_k, t_undi,q_undi
2278	o ,wake_omgbdth,wake_dp_omgb
2279	o ,wake_dtKE,wake_dqKE
2280	o ,wake_dtPBL,wake_dqPBL
2281	o ,wake_omg,wake_dp_deltomg
2282	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2283	o ,wake_ddeltat,wake_ddeltaq)
2284	c
2285	!-----------------------------------------------------------------------------------------
2286	! ajout des tendances des poches froides
2287	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2288	! coherent avec les autres d_t_...
2289	d_t_wake(:,:)=dt_wake(:,:)*dtime
2290	d_q_wake(:,:)=dq_wake(:,:)*dtime
2291	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2292	!-----------------------------------------------------------------------------------------
2293
2294	endif
2295	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2296	c
2297	c===================================================================
2298	c Convection seche (thermiques ou ajustement)
2299	c===================================================================
2300	c
2301	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2302	s ,seuil_inversion,weak_inversion,dthmin)
2303
2304
2305
2306	d_t_ajsb(:,:)=0.
2307	d_q_ajsb(:,:)=0.
2308	d_t_ajs(:,:)=0.
2309	d_u_ajs(:,:)=0.
2310	d_v_ajs(:,:)=0.
2311	d_q_ajs(:,:)=0.
2312	clwcon0th(:,:)=0.
2313	c
2314	fm_therm(:,:)=0.
2315	entr_therm(:,:)=0.
2316	detr_therm(:,:)=0.
2317	c
2318	IF(prt_level>9)WRITE(lunout,*)
2319	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2320	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2321	if(iflag_thermals.lt.0) then
2322	c Rien
2323	c ====
2324	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2325
2326
2327	else
2328
2329	c Thermiques
2330	c ==========
2331	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2332	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2333
2334
2335	if (iflag_thermals.gt.1) then
2336	call calltherm(pdtphys
2337	s ,pplay,paprs,pphi,weak_inversion
2338	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2339	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2340	s ,fm_therm,entr_therm,detr_therm
2341	s ,zqasc,clwcon0th,lmax_th,ratqscth
2342	s ,ratqsdiff,zqsatth
2343	con rajoute ale et alp, et les caracteristiques de la couche alim
2344	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
2345	endif
2346
2347
2348	c Ajustement sec
2349	c ==============
2350
2351	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2352	! a partir du sommet des thermiques.
2353	! Dans le cas contraire, on demarre au niveau 1.
2354
2355	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2356
2357	if(iflag_thermals.eq.0) then
2358	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2359	limbas(:)=1
2360	else
2361	limbas(:)=lmax_th(:)
2362	endif
2363
2364	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2365	! pour des test de convergence numerique.
2366	! Le nouveau ajsec est a priori mieux, meme pour le cas
2367	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2368	! non nulles numeriquement pour des mailles non concernees.
2369
2370	if (iflag_thermals.eq.0) then
2371	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2372	s , d_t_ajsb, d_q_ajsb)
2373	else
2374	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2375	s , d_t_ajsb, d_q_ajsb)
2376	endif
2377
2378	!-----------------------------------------------------------------------------------------
2379	! ajout des tendances de l'ajustement sec ou des thermiques
2380	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2381	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2382	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2383
2384	!-----------------------------------------------------------------------------------------
2385
2386	endif
2387
2388	endif
2389	c
2390	c===================================================================
2391	cIM
2392	IF (ip_ebil_phy.ge.2) THEN
2393	ztit='after dry_adjust'
2394	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2395	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2396	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2397	END IF
2398
2399
2400	c-------------------------------------------------------------------------
2401	c Caclul des ratqs
2402	c-------------------------------------------------------------------------
2403
2404	c print*,'calcul des ratqs'
2405	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2406	c ----------------
2407	c on ecrase le tableau ratqsc calcule par clouds_gno
2408	if (iflag_cldcon.eq.1) then
2409	do k=1,klev
2410	do i=1,klon
2411	if(ptconv(i,k)) then
2412	ratqsc(i,k)=ratqsbas
2413	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2414	else
2415	ratqsc(i,k)=0.
2416	endif
2417	enddo
2418	enddo
2419
2420	c-----------------------------------------------------------------------
2421	c par nversion de la fonction log normale
2422	c-----------------------------------------------------------------------
2423	else if (iflag_cldcon.eq.4) then
2424	ptconvth(:,:)=.false.
2425	ratqsc(:,:)=0.
2426	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2427	call clouds_gno
2428	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2429	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2430
2431	c-----------------------------------------------------------------------
2432	c par calcul direct de l'ecart-type
2433	c-----------------------------------------------------------------------
2434
2435	else if (iflag_cldcon>=5) then
2436	wmax_th(:)=0.
2437	zmax_th(:)=0.
2438	do k=1,klev
2439	do i=1,klon
2440	wmax_th(i)=max(wmax_th(i),zw2(i,k))
2441	if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
2442	enddo
2443	enddo
2444	tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
2445	print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
2446
2447	c On impose que l'air autour de la fraction couverte par le thermique
2448	c plus son air detraine durant tau_overturning_th soit superieur
2449	c a 0.1 q_seri
2450	zz=0.1
2451	do k=1,klev
2452	do i=1,klon
2453	lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
2454	s tau_overturning_th(i)*detr_therm(i,k)
2455	s *rg/(paprs(i,k)-paprs(i,k+1))
2456	znum=(1.-zz)*q_seri(i,k)
2457	zden=zqasc(i,k)-zz*q_seri(i,k)
2458	if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
2459	lambda_th(i,k)=min(lambda_th(i,k),0.9)
2460	enddo
2461	enddo
2462
2463	if(iflag_cldcon==5) then
2464	do k=1,klev
2465	do i=1,klon
2466	ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
2467	s abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
2468	enddo
2469	enddo
2470	else if(iflag_cldcon==6) then
2471	do k=1,klev
2472	do i=1,klon
2473	ratqsc(i,k)=sqrt(lambda_th(i,k))*
2474	s (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
2475	enddo
2476	enddo
2477	endif
2478
2479	endif
2480
2481	c ratqs stables
2482	c -------------
2483
2484	if (iflag_ratqs.eq.0) then
2485
2486	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2487	do k=1,klev
2488	do i=1, klon
2489	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2490	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2491	enddo
2492	enddo
2493
2494	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2495	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2496	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2497	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2498	! Il s'agit de differents tests dans la phase de reglage du modele
2499	! avec thermiques.
2500
2501	else if (iflag_ratqs.eq.1) then
2502
2503	do k=1,klev
2504	do i=1, klon
2505	if (pplay(i,k).ge.60000.) then
2506	ratqss(i,k)=ratqsbas
2507	else if ((pplay(i,k).ge.30000.).and.
2508	s (pplay(i,k).lt.60000.)) then
2509	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2510	s (60000.-pplay(i,k))/(60000.-30000.)
2511	else
2512	ratqss(i,k)=ratqshaut
2513	endif
2514	enddo
2515	enddo
2516
2517	else
2518
2519	do k=1,klev
2520	do i=1, klon
2521	if (pplay(i,k).ge.60000.) then
2522	ratqss(i,k)=ratqsbas
2523	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2524	else if ((pplay(i,k).ge.30000.).and.
2525	s (pplay(i,k).lt.60000.)) then
2526	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2527	s (60000.-pplay(i,k))/(60000.-30000.)
2528	else
2529	ratqss(i,k)=ratqshaut
2530	endif
2531	enddo
2532	enddo
2533	endif
2534
2535
2536
2537
2538	c ratqs final
2539	c -----------
2540
2541	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2542	s .or.iflag_cldcon.ge.4) then
2543
2544	! On ajoute une constante au ratqsc*2 pour tenir compte de
2545	! fluctuations turbulentes de petite echelle
2546
2547	do k=1,klev
2548	do i=1,klon
2549	if ((fm_therm(i,k).gt.1.e-10)) then
2550	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2551	endif
2552	enddo
2553	enddo
2554
2555	! les ratqs sont une combinaison de ratqss et ratqsc
2556	! print*,'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
2557
2558	if (tau_ratqs>1.e-10) then
2559	facteur=exp(-pdtphys/tau_ratqs)
2560	else
2561	facteur=0.
2562	endif
2563	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2564	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2565	! FH 22/09/2009
2566	! La ligne ci-dessous faisait osciller le modele et donnait une solution
2567	! assymptotique bidon et dÃ©pendant fortement du pas de temps.
2568	! ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2569	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2570	ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
2571	else
2572	! on ne prend que le ratqs stable pour fisrtilp
2573	ratqs(:,:)=ratqss(:,:)
2574	endif
2575
2576
2577	c
2578	c Appeler le processus de condensation a grande echelle
2579	c et le processus de precipitation
2580	c-------------------------------------------------------------------------
2581	CALL fisrtilp(dtime,paprs,pplay,
2582	. t_seri, q_seri,ptconv,ratqs,
2583	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2584	. rain_lsc, snow_lsc,
2585	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2586	. frac_impa, frac_nucl,
2587	. prfl, psfl, rhcl)
2588
2589	WHERE (rain_lsc < 0) rain_lsc = 0.
2590	WHERE (snow_lsc < 0) snow_lsc = 0.
2591	!-----------------------------------------------------------------------------------------
2592	! ajout des tendances de la diffusion turbulente
2593	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2594	!-----------------------------------------------------------------------------------------
2595	DO k = 1, klev
2596	DO i = 1, klon
2597	cldfra(i,k) = rneb(i,k)
2598	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2599	ENDDO
2600	ENDDO
2601	IF (check) THEN
2602	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2603	WRITE(lunout,*)"apresilp=", za
2604	zx_t = 0.0
2605	za = 0.0
2606	DO i = 1, klon
2607	za = za + airephy(i)/FLOAT(klon)
2608	zx_t = zx_t + (rain_lsc(i)
2609	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2610	ENDDO
2611	zx_t = zx_t/za*dtime
2612	WRITE(lunout,*)"Precip=", zx_t
2613	ENDIF
2614	cIM
2615	IF (ip_ebil_phy.ge.2) THEN
2616	ztit='after fisrt'
2617	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2618	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2619	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2620	call diagphy(airephy,ztit,ip_ebil_phy
2621	e , zero_v, zero_v, zero_v, zero_v, zero_v
2622	e , zero_v, rain_lsc, snow_lsc, ztsol
2623	e , d_h_vcol, d_qt, d_ec
2624	s , fs_bound, fq_bound )
2625	END IF
2626
2627	if (mydebug) then
2628	call writefield_phy('u_seri',u_seri,llm)
2629	call writefield_phy('v_seri',v_seri,llm)
2630	call writefield_phy('t_seri',t_seri,llm)
2631	call writefield_phy('q_seri',q_seri,llm)
2632	endif
2633
2634	c
2635	c-------------------------------------------------------------------
2636	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2637	c-------------------------------------------------------------------
2638
2639	c 1. NUAGES CONVECTIFS
2640	c
2641	cIM cf FH
2642	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2643	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2644	snow_tiedtke=0.
2645	c print*,'avant calcul de la pseudo precip '
2646	c print*,'iflag_cldcon',iflag_cldcon
2647	if (iflag_cldcon.eq.-1) then
2648	rain_tiedtke=rain_con
2649	else
2650	c print*,'calcul de la pseudo precip '
2651	rain_tiedtke=0.
2652	c print*,'calcul de la pseudo precip 0'
2653	do k=1,klev
2654	do i=1,klon
2655	if (d_q_con(i,k).lt.0.) then
2656	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2657	s *(paprs(i,k)-paprs(i,k+1))/rg
2658	endif
2659	enddo
2660	enddo
2661	endif
2662	c
2663	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2664	c
2665
2666	c Nuages diagnostiques pour Tiedtke
2667	CALL diagcld1(paprs,pplay,
2668	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2669	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2670	. diafra,dialiq)
2671	DO k = 1, klev
2672	DO i = 1, klon
2673	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2674	cldliq(i,k) = dialiq(i,k)
2675	cldfra(i,k) = diafra(i,k)
2676	ENDIF
2677	ENDDO
2678	ENDDO
2679
2680	ELSE IF (iflag_cldcon.ge.3) THEN
2681	c On prend pour les nuages convectifs le max du calcul de la
2682	c convection et du calcul du pas de temps precedent diminue d'un facteur
2683	c facttemps
2684	facteur = pdtphys *facttemps
2685	do k=1,klev
2686	do i=1,klon
2687	rnebcon(i,k)=rnebcon(i,k)*facteur
2688	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2689	s then
2690	rnebcon(i,k)=rnebcon0(i,k)
2691	clwcon(i,k)=clwcon0(i,k)
2692	endif
2693	enddo
2694	enddo
2695
2696	c
2697	cjq - introduce the aerosol direct and first indirect radiative forcings
2698	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2699	IF (ok_ade.OR.ok_aie) THEN
2700	IF (.NOT. aerosol_couple)
2701	& CALL readaerosol_optic(
2702	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
2703	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
2704	& mass_solu_aero, mass_solu_aero_pi,
2705	& tau_aero, piz_aero, cg_aero,
2706	& tausum_aero, tau3d_aero)
2707	ELSE
2708	tau_aero(:,:,:,:) = 0.
2709	piz_aero(:,:,:,:) = 0.
2710	cg_aero(:,:,:,:) = 0.
2711	ENDIF
2712
2713	cIM calcul nuages par le simulateur ISCCP
2714	c
2715	#ifdef histISCCP
2716	IF (ok_isccp) THEN
2717	c
2718	cIM lecture invtau, tautab des fichiers formattes
2719	c
2720	IF (debut) THEN
2721	c$OMP MASTER
2722	c
2723	open(99,file='tautab.formatted', FORM='FORMATTED')
2724	read(99,'(f30.20)') tautab_omp
2725	close(99)
2726	c
2727	open(99,file='invtau.formatted',form='FORMATTED')
2728	read(99,'(i10)') invtau_omp
2729
2730	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
2731	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2732
2733	close(99)
2734	c$OMP END MASTER
2735	c$OMP BARRIER
2736	tautab=tautab_omp
2737	invtau=invtau_omp
2738	c
2739	ENDIF !debut
2740	c
2741	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2742	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2743	#include "calcul_simulISCCP.h"
2744	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2745	ENDIF !ok_isccp
2746	#endif
2747
2748	c On prend la somme des fractions nuageuses et des contenus en eau
2749	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2750	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2751
2752	ENDIF
2753	c
2754	c 2. NUAGES STARTIFORMES
2755	c
2756	IF (ok_stratus) THEN
2757	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2758	DO k = 1, klev
2759	DO i = 1, klon
2760	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2761	cldliq(i,k) = dialiq(i,k)
2762	cldfra(i,k) = diafra(i,k)
2763	ENDIF
2764	ENDDO
2765	ENDDO
2766	ENDIF
2767	c
2768	c Precipitation totale
2769	c
2770	DO i = 1, klon
2771	rain_fall(i) = rain_con(i) + rain_lsc(i)
2772	snow_fall(i) = snow_con(i) + snow_lsc(i)
2773	ENDDO
2774	cIM
2775	IF (ip_ebil_phy.ge.2) THEN
2776	ztit="after diagcld"
2777	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2778	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2779	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2780	END IF
2781	c
2782	c Calculer l'humidite relative pour diagnostique
2783	c
2784	DO k = 1, klev
2785	DO i = 1, klon
2786	zx_t = t_seri(i,k)
2787	IF (thermcep) THEN
2788	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2789	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2790	zx_qs = MIN(0.5,zx_qs)
2791	zcor = 1./(1.-retv*zx_qs)
2792	zx_qs = zx_qs*zcor
2793	ELSE
2794	IF (zx_t.LT.t_coup) THEN
2795	zx_qs = qsats(zx_t)/pplay(i,k)
2796	ELSE
2797	zx_qs = qsatl(zx_t)/pplay(i,k)
2798	ENDIF
2799	ENDIF
2800	zx_rh(i,k) = q_seri(i,k)/zx_qs
2801	zqsat(i,k)=zx_qs
2802	ENDDO
2803	ENDDO
2804
2805	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2806	c equivalente a 2m (tpote) pour diagnostique
2807	c
2808	DO i = 1, klon
2809	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2810	IF (thermcep) THEN
2811	IF(zt2m(i).LT.RTT) then
2812	Lheat=RLSTT
2813	ELSE
2814	Lheat=RLVTT
2815	ENDIF
2816	ELSE
2817	IF (zt2m(i).LT.RTT) THEN
2818	Lheat=RLSTT
2819	ELSE
2820	Lheat=RLVTT
2821	ENDIF
2822	ENDIF
2823	tpote(i) = tpot(i)*
2824	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2825	ENDDO
2826
2827	IF (config_inca /= 'none') THEN
2828	#ifdef INCA
2829	CALL VTe(VTphysiq)
2830	CALL VTb(VTinca)
2831	calday = FLOAT(days_elapsed + 1) + jH_cur
2832
2833	call chemtime(itap+itau_phy-1, date0, dtime)
2834	IF (config_inca == 'aero') THEN
2835	CALL AEROSOL_METEO_CALC(
2836	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2837	$ prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
2838	END IF
2839
2840	zxsnow_dummy(:) = 0.0
2841
2842	CALL chemhook_begin (calday,
2843	$ days_elapsed+1,
2844	$ jH_cur,
2845	$ pctsrf(1,1),
2846	$ rlat,
2847	$ rlon,
2848	$ airephy,
2849	$ paprs,
2850	$ pplay,
2851	$ coefh,
2852	$ pphi,
2853	$ t_seri,
2854	$ u,
2855	$ v,
2856	$ wo(:, :, 1),
2857	$ q_seri,
2858	$ zxtsol,
2859	$ zxsnow_dummy,
2860	$ solsw,
2861	$ albsol1,
2862	$ rain_fall,
2863	$ snow_fall,
2864	$ itop_con,
2865	$ ibas_con,
2866	$ cldfra,
2867	$ iim,
2868	$ jjm,
2869	$ tr_seri,
2870	$ ftsol,
2871	$ paprs,
2872	$ cdragh,
2873	$ cdragm,
2874	$ pctsrf,
2875	$ pdtphys,
2876	$ itap)
2877
2878	CALL VTe(VTinca)
2879	CALL VTb(VTphysiq)
2880	#endif
2881	END IF !config_inca /= 'none'
2882	c
2883	c Calculer les parametres optiques des nuages et quelques
2884	c parametres pour diagnostiques:
2885	c
2886
2887	IF (aerosol_couple) THEN
2888	mass_solu_aero(:,:) = ccm(:,:,1)
2889	mass_solu_aero_pi(:,:) = ccm(:,:,2)
2890	END IF
2891
2892	if (ok_newmicro) then
2893	CALL newmicro (paprs, pplay,ok_newmicro,
2894	. t_seri, cldliq, cldfra, cldtau, cldemi,
2895	. cldh, cldl, cldm, cldt, cldq,
2896	. flwp, fiwp, flwc, fiwc,
2897	e ok_aie,
2898	e mass_solu_aero, mass_solu_aero_pi,
2899	e bl95_b0, bl95_b1,
2900	s cldtaupi, re, fl, ref_liq, ref_ice)
2901	else
2902	CALL nuage (paprs, pplay,
2903	. t_seri, cldliq, cldfra, cldtau, cldemi,
2904	. cldh, cldl, cldm, cldt, cldq,
2905	e ok_aie,
2906	e mass_solu_aero, mass_solu_aero_pi,
2907	e bl95_b0, bl95_b1,
2908	s cldtaupi, re, fl)
2909
2910	endif
2911	c
2912	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2913	c
2914	IF (MOD(itaprad,radpas).EQ.0) THEN
2915
2916	DO i = 1, klon
2917	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2918	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2919	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2920	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2921	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2922	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2923	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2924	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2925	ENDDO
2926
2927	if (mydebug) then
2928	call writefield_phy('u_seri',u_seri,llm)
2929	call writefield_phy('v_seri',v_seri,llm)
2930	call writefield_phy('t_seri',t_seri,llm)
2931	call writefield_phy('q_seri',q_seri,llm)
2932	endif
2933
2934	IF (aerosol_couple) THEN
2935	#ifdef INCA
2936	CALL radlwsw_inca
2937	e (kdlon,kflev,dist, rmu0, fract, solaire,
2938	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2939	e wo(:, :, 1),
2940	e cldfra, cldemi, cldtau,
2941	s heat,heat0,cool,cool0,radsol,albpla,
2942	s topsw,toplw,solsw,sollw,
2943	s sollwdown,
2944	s topsw0,toplw0,solsw0,sollw0,
2945	s lwdn0, lwdn, lwup0, lwup,
2946	s swdn0, swdn, swup0, swup,
2947	e ok_ade, ok_aie,
2948	e tau_aero, piz_aero, cg_aero,
2949	s topswad_aero, solswad_aero,
2950	s topswad0_aero, solswad0_aero,
2951	s topsw_aero, topsw0_aero,
2952	s solsw_aero, solsw0_aero,
2953	e cldtaupi,
2954	s topswai_aero, solswai_aero)
2955
2956	#endif
2957	ELSE
2958
2959	CALL radlwsw
2960	e (dist, rmu0, fract,
2961	e paprs, pplay,zxtsol,albsol1, albsol2,
2962	e t_seri,q_seri,wo,
2963	e cldfra, cldemi, cldtau,
2964	e ok_ade, ok_aie,
2965	e tau_aero, piz_aero, cg_aero,
2966	e cldtaupi,new_aod,
2967	e zqsat, flwc, fiwc,
2968	s heat,heat0,cool,cool0,radsol,albpla,
2969	s topsw,toplw,solsw,sollw,
2970	s sollwdown,
2971	s topsw0,toplw0,solsw0,sollw0,
2972	s lwdn0, lwdn, lwup0, lwup,
2973	s swdn0, swdn, swup0, swup,
2974	s topswad_aero, solswad_aero,
2975	s topswai_aero, solswai_aero,
2976	o topswad0_aero, solswad0_aero,
2977	o topsw_aero, topsw0_aero,
2978	o solsw_aero, solsw0_aero,
2979	o topswcf_aero, solswcf_aero)
2980
2981
2982	ENDIF ! aerosol_couple
2983	itaprad = 0
2984	ENDIF ! MOD(itaprad,radpas)
2985	itaprad = itaprad + 1
2986
2987	if (iflag_radia.eq.0) then
2988	print *,'--------------------------------------------------'
2989	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
2990	print *,'>>>> heat et cool mis a zero '
2991	print *,'--------------------------------------------------'
2992	heat=0.
2993	cool=0.
2994	endif
2995
2996	c
2997	c Ajouter la tendance des rayonnements (tous les pas)
2998	c
2999	DO k = 1, klev
3000	DO i = 1, klon
3001	t_seri(i,k) = t_seri(i,k)
3002	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
3003	ENDDO
3004	ENDDO
3005	c
3006	if (mydebug) then
3007	call writefield_phy('u_seri',u_seri,llm)
3008	call writefield_phy('v_seri',v_seri,llm)
3009	call writefield_phy('t_seri',t_seri,llm)
3010	call writefield_phy('q_seri',q_seri,llm)
3011	endif
3012
3013	cIM
3014	IF (ip_ebil_phy.ge.2) THEN
3015	ztit='after rad'
3016	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3017	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3018	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3019	call diagphy(airephy,ztit,ip_ebil_phy
3020	e , topsw, toplw, solsw, sollw, zero_v
3021	e , zero_v, zero_v, zero_v, ztsol
3022	e , d_h_vcol, d_qt, d_ec
3023	s , fs_bound, fq_bound )
3024	END IF
3025	c
3026	c
3027	c Calculer l'hydrologie de la surface
3028	c
3029	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3030	c . agesno, ftsol,fqsurf,fsnow, ruis)
3031	c
3032
3033	c
3034	c Calculer le bilan du sol et la derive de temperature (couplage)
3035	c
3036	DO i = 1, klon
3037	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3038	c a la demande de JLD
3039	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3040	ENDDO
3041	c
3042	cmoddeblott(jan95)
3043	c Appeler le programme de parametrisation de l'orographie
3044	c a l'echelle sous-maille:
3045	c
3046	IF (ok_orodr) THEN
3047	c
3048	c selection des points pour lesquels le shema est actif:
3049	igwd=0
3050	DO i=1,klon
3051	itest(i)=0
3052	c IF ((zstd(i).gt.10.0)) THEN
3053	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3054	itest(i)=1
3055	igwd=igwd+1
3056	idx(igwd)=i
3057	ENDIF
3058	ENDDO
3059	c igwdim=MAX(1,igwd)
3060	c
3061	IF (ok_strato) THEN
3062
3063	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3064	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3065	e igwd,idx,itest,
3066	e t_seri, u_seri, v_seri,
3067	s zulow, zvlow, zustrdr, zvstrdr,
3068	s d_t_oro, d_u_oro, d_v_oro)
3069
3070	ELSE
3071	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3072	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3073	e igwd,idx,itest,
3074	e t_seri, u_seri, v_seri,
3075	s zulow, zvlow, zustrdr, zvstrdr,
3076	s d_t_oro, d_u_oro, d_v_oro)
3077	ENDIF
3078	c
3079	c ajout des tendances
3080	!-----------------------------------------------------------------------------------------
3081	! ajout des tendances de la trainee de l'orographie
3082	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3083	!-----------------------------------------------------------------------------------------
3084	c
3085	ENDIF ! fin de test sur ok_orodr
3086	c
3087	if (mydebug) then
3088	call writefield_phy('u_seri',u_seri,llm)
3089	call writefield_phy('v_seri',v_seri,llm)
3090	call writefield_phy('t_seri',t_seri,llm)
3091	call writefield_phy('q_seri',q_seri,llm)
3092	endif
3093
3094	IF (ok_orolf) THEN
3095	c
3096	c selection des points pour lesquels le shema est actif:
3097	igwd=0
3098	DO i=1,klon
3099	itest(i)=0
3100	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3101	itest(i)=1
3102	igwd=igwd+1
3103	idx(igwd)=i
3104	ENDIF
3105	ENDDO
3106	c igwdim=MAX(1,igwd)
3107	c
3108	IF (ok_strato) THEN
3109
3110	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3111	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3112	e igwd,idx,itest,
3113	e t_seri, u_seri, v_seri,
3114	s zulow, zvlow, zustrli, zvstrli,
3115	s d_t_lif, d_u_lif, d_v_lif )
3116
3117	ELSE
3118	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3119	e rlat,zmea,zstd,zpic,
3120	e itest,
3121	e t_seri, u_seri, v_seri,
3122	s zulow, zvlow, zustrli, zvstrli,
3123	s d_t_lif, d_u_lif, d_v_lif)
3124	ENDIF
3125	c
3126	!-----------------------------------------------------------------------------------------
3127	! ajout des tendances de la portance de l'orographie
3128	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3129	!-----------------------------------------------------------------------------------------
3130	c
3131	ENDIF ! fin de test sur ok_orolf
3132	C HINES GWD PARAMETRIZATION
3133
3134	IF (ok_hines) then
3135
3136	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3137	i rlat,t_seri,u_seri,v_seri,
3138	o zustrhi,zvstrhi,
3139	o d_t_hin, d_u_hin, d_v_hin)
3140	c
3141	c ajout des tendances
3142	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3143
3144	ENDIF
3145	c
3146
3147	c
3148	cIM cf. FLott BEG
3149	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3150
3151	if (mydebug) then
3152	call writefield_phy('u_seri',u_seri,llm)
3153	call writefield_phy('v_seri',v_seri,llm)
3154	call writefield_phy('t_seri',t_seri,llm)
3155	call writefield_phy('q_seri',q_seri,llm)
3156	endif
3157
3158	DO i = 1, klon
3159	zustrph(i)=0.
3160	zvstrph(i)=0.
3161	ENDDO
3162	DO k = 1, klev
3163	DO i = 1, klon
3164	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3165	c (paprs(i,k)-paprs(i,k+1))/rg
3166	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3167	c (paprs(i,k)-paprs(i,k+1))/rg
3168	ENDDO
3169	ENDDO
3170	c
3171	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3172	c
3173	IF (is_sequential) THEN
3174
3175	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3176	C ra,rg,romega,
3177	C rlat,rlon,pphis,
3178	C zustrdr,zustrli,zustrph,
3179	C zvstrdr,zvstrli,zvstrph,
3180	C paprs,u,v,
3181	C aam, torsfc)
3182	ENDIF
3183	cIM cf. FLott END
3184	cIM
3185	IF (ip_ebil_phy.ge.2) THEN
3186	ztit='after orography'
3187	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3188	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3189	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3190	END IF
3191	c
3192	c
3193	!====================================================================
3194	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3195	!====================================================================
3196	! Abderrahmane 24.08.09
3197
3198	IF (ok_cosp) THEN
3199	! adeclarer
3200	#ifdef CPP_COSP
3201	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3202
3203	print*,'freq_cosp',freq_cosp
3204	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3205	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3206	! s ref_liq,ref_ice
3207	call phys_cosp(itap,dtime,freq_cosp,
3208	$ ecrit_mth,ecrit_day,ecrit_hf,overlap,
3209	$ klon,klev,rlon,rlat,presnivs,
3210	$ ref_liq,ref_ice,
3211	$ pctsrf(:,is_ter)+pctsrf(:,is_lic),
3212	$ zu10m,zv10m,
3213	$ zphi,paprs(:,1:klev),pplay,zxtsol,t_seri,
3214	$ qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc,
3215	$ prfl(:,1:klev),psfl(:,1:klev),
3216	$ pmflxr(:,1:klev),pmflxs(:,1:klev),
3217	$ mr_ozone,cldtau, cldemi)
3218
3219	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3220	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3221	! M clMISR,
3222	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3223	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3224
3225	ENDIF
3226
3227	#endif
3228	ENDIF !ok_cosp
3229	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3230	cAA
3231	cAA Installation de l'interface online-offline pour traceurs
3232	cAA
3233	c====================================================================
3234	c Calcul des tendances traceurs
3235	c====================================================================
3236	C
3237
3238	call phytrac (
3239	I itap, days_elapsed+1, jH_cur, debut,
3240	I lafin, dtime, u, v, t,
3241	I paprs, pplay, pmfu, pmfd,
3242	I pen_u, pde_u, pen_d, pde_d,
3243	I cdragh, coefh, fm_therm, entr_therm,
3244	I u1, v1, ftsol, pctsrf,
3245	I rlat, frac_impa, frac_nucl,rlon,
3246	I presnivs, pphis, pphi, albsol1,
3247	I qx(:,:,ivap),rhcl, cldfra, rneb,
3248	I diafra, cldliq, itop_con, ibas_con,
3249	I pmflxr, pmflxs, prfl, psfl,
3250	I da, phi, mp, upwd,
3251	I dnwd, aerosol_couple, flxmass_w,
3252	I tau_aero, piz_aero, cg_aero, ccm,
3253	I rfname,
3254	O tr_seri)
3255
3256	IF (offline) THEN
3257
3258	print*,'Attention on met a 0 les thermiques pour phystoke'
3259	call phystokenc (
3260	I nlon,klev,pdtphys,rlon,rlat,
3261	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3262	I fm_therm,entr_therm,
3263	I cdragh,coefh,u1,v1,ftsol,pctsrf,
3264	I frac_impa, frac_nucl,
3265	I pphis,airephy,dtime,itap)
3266
3267
3268	ENDIF
3269
3270	c
3271	c Calculer le transport de l'eau et de l'energie (diagnostique)
3272	c
3273	CALL transp (paprs,zxtsol,
3274	e t_seri, q_seri, u_seri, v_seri, zphi,
3275	s ve, vq, ue, uq)
3276	c
3277	cIM global posePB BEG
3278	IF(1.EQ.0) THEN
3279	c
3280	CALL transp_lay (paprs,zxtsol,
3281	e t_seri, q_seri, u_seri, v_seri, zphi,
3282	s ve_lay, vq_lay, ue_lay, uq_lay)
3283	c
3284	ENDIF !(1.EQ.0) THEN
3285	cIM global posePB END
3286	c Accumuler les variables a stocker dans les fichiers histoire:
3287	c
3288	c+jld ec_conser
3289	DO k = 1, klev
3290	DO i = 1, klon
3291	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3292	d_t_ec(i,k)=0.5/ZRCPD
3293	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3294	ENDDO
3295	ENDDO
3296
3297	DO k = 1, klev
3298	DO i = 1, klon
3299	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3300	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3301	END DO
3302	END DO
3303	c-jld ec_conser
3304	cIM
3305	IF (ip_ebil_phy.ge.1) THEN
3306	ztit='after physic'
3307	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3308	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3309	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3310	C Comme les tendances de la physique sont ajoute dans la dynamique,
3311	C on devrait avoir que la variation d'entalpie par la dynamique
3312	C est egale a la variation de la physique au pas de temps precedent.
3313	C Donc la somme de ces 2 variations devrait etre nulle.
3314
3315	call diagphy(airephy,ztit,ip_ebil_phy
3316	e , topsw, toplw, solsw, sollw, sens
3317	e , evap, rain_fall, snow_fall, ztsol
3318	e , d_h_vcol, d_qt, d_ec
3319	s , fs_bound, fq_bound )
3320	C
3321	d_h_vcol_phy=d_h_vcol
3322	C
3323	END IF
3324	C
3325	c=======================================================================
3326	c SORTIES
3327	c=======================================================================
3328
3329	cIM Interpolation sur les niveaux de pression du NMC
3330	c -------------------------------------------------
3331	c
3332	#include "calcul_STDlev.h"
3333	twriteSTD(:,:,1)=tsumSTD(:,:,2)
3334	qwriteSTD(:,:,1)=qsumSTD(:,:,2)
3335	rhwriteSTD(:,:,1)=rhsumSTD(:,:,2)
3336	phiwriteSTD(:,:,1)=phisumSTD(:,:,2)
3337	uwriteSTD(:,:,1)=usumSTD(:,:,2)
3338	vwriteSTD(:,:,1)=vsumSTD(:,:,2)
3339	wwriteSTD(:,:,1)=wsumSTD(:,:,2)
3340
3341	twriteSTD(:,:,2)=tsumSTD(:,:,1)
3342	qwriteSTD(:,:,2)=qsumSTD(:,:,1)
3343	rhwriteSTD(:,:,2)=rhsumSTD(:,:,1)
3344	phiwriteSTD(:,:,2)=phisumSTD(:,:,1)
3345	uwriteSTD(:,:,2)=usumSTD(:,:,1)
3346	vwriteSTD(:,:,2)=vsumSTD(:,:,1)
3347	wwriteSTD(:,:,2)=wsumSTD(:,:,1)
3348
3349	twriteSTD(:,:,3)=tlevSTD(:,:)
3350	qwriteSTD(:,:,3)=qlevSTD(:,:)
3351	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3352	phiwriteSTD(:,:,3)=philevSTD(:,:)
3353	uwriteSTD(:,:,3)=ulevSTD(:,:)
3354	vwriteSTD(:,:,3)=vlevSTD(:,:)
3355	wwriteSTD(:,:,3)=wlevSTD(:,:)
3356
3357	twriteSTD(:,:,4)=tlevSTD(:,:)
3358	qwriteSTD(:,:,4)=qlevSTD(:,:)
3359	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3360	phiwriteSTD(:,:,4)=philevSTD(:,:)
3361	uwriteSTD(:,:,4)=ulevSTD(:,:)
3362	vwriteSTD(:,:,4)=vlevSTD(:,:)
3363	wwriteSTD(:,:,4)=wlevSTD(:,:)
3364	c
3365	c slp sea level pressure
3366	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3367	c
3368	ccc prw = eau precipitable
3369	DO i = 1, klon
3370	prw(i) = 0.
3371	DO k = 1, klev
3372	prw(i) = prw(i) +
3373	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3374	ENDDO
3375	ENDDO
3376	c
3377	cIM initialisation + calculs divers diag AMIP2
3378	c
3379	#include "calcul_divers.h"
3380	c
3381	IF (config_inca /= 'none') THEN
3382	#ifdef INCA
3383	CALL VTe(VTphysiq)
3384	CALL VTb(VTinca)
3385
3386	CALL chemhook_end (
3387	$ dtime,
3388	$ pplay,
3389	$ t_seri,
3390	$ tr_seri,
3391	$ nbtr,
3392	$ paprs,
3393	$ q_seri,
3394	$ airephy,
3395	$ pphi,
3396	$ pphis,
3397	$ zx_rh)
3398
3399	CALL VTe(VTinca)
3400	CALL VTb(VTphysiq)
3401	#endif
3402	END IF
3403
3404	c=============================================================
3405	c
3406	c Convertir les incrementations en tendances
3407	c
3408	if (mydebug) then
3409	call writefield_phy('u_seri',u_seri,llm)
3410	call writefield_phy('v_seri',v_seri,llm)
3411	call writefield_phy('t_seri',t_seri,llm)
3412	call writefield_phy('q_seri',q_seri,llm)
3413	endif
3414
3415	DO k = 1, klev
3416	DO i = 1, klon
3417	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3418	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3419	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3420	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3421	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3422	ENDDO
3423	ENDDO
3424	c
3425	IF (nqtot.GE.3) THEN
3426	DO iq = 3, nqtot
3427	DO k = 1, klev
3428	DO i = 1, klon
3429	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3430	ENDDO
3431	ENDDO
3432	ENDDO
3433	ENDIF
3434	c
3435	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3436	cIM global posePB#include "write_bilKP_ins.h"
3437	cIM global posePB#include "write_bilKP_ave.h"
3438	c
3439
3440	c Sauvegarder les valeurs de t et q a la fin de la physique:
3441	c
3442	DO k = 1, klev
3443	DO i = 1, klon
3444	u_ancien(i,k) = u_seri(i,k)
3445	v_ancien(i,k) = v_seri(i,k)
3446	t_ancien(i,k) = t_seri(i,k)
3447	q_ancien(i,k) = q_seri(i,k)
3448	ENDDO
3449	ENDDO
3450	c
3451	!==========================================================================
3452	! Sorties des tendances pour un point particulier
3453	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3454	! pour le debug
3455	! La valeur de igout est attribuee plus haut dans le programme
3456	!==========================================================================
3457
3458	if (prt_level.ge.1) then
3459	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3460	write(lunout,*)
3461	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3462	write(lunout,*)
3463	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
3464	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3465	s pctsrf(igout,is_sic)
3466	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3467	do k=1,klev
3468	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3469	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3470	s d_t_eva(igout,k)
3471	enddo
3472	write(lunout,*) 'cool,heat'
3473	do k=1,klev
3474	write(lunout,*) cool(igout,k),heat(igout,k)
3475	enddo
3476
3477	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3478	do k=1,klev
3479	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3480	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3481	enddo
3482
3483	write(lunout,*) 'd_ps ',d_ps(igout)
3484	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3485	do k=1,klev
3486	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3487	s d_qx(igout,k,1),d_qx(igout,k,2)
3488	enddo
3489	endif
3490
3491	!==========================================================================
3492
3493	c============================================================
3494	c Calcul de la temperature potentielle
3495	c============================================================
3496	DO k = 1, klev
3497	DO i = 1, klon
3498	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3499	ENDDO
3500	ENDDO
3501	c
3502
3503	c 22.03.04 BEG
3504	c=============================================================
3505	c Ecriture des sorties
3506	c=============================================================
3507	#ifdef CPP_IOIPSL
3508
3509	c Recupere des varibles calcule dans differents modules
3510	c pour ecriture dans histxxx.nc
3511
3512	! Get some variables from module fonte_neige_mod
3513	CALL fonte_neige_get_vars(pctsrf,
3514	. zxfqcalving, zxfqfonte, zxffonte)
3515
3516
3517	#include "phys_output_write.h"
3518
3519	#ifdef histISCCP
3520	#include "write_histISCCP.h"
3521	#endif
3522
3523	#ifdef histmthNMC
3524	#include "write_histmthNMC.h"
3525	#endif
3526
3527	#include "write_histday_seri.h"
3528
3529	#include "write_paramLMDZ_phy.h"
3530
3531	#endif
3532
3533	c 22.03.04 END
3534	c
3535	c====================================================================
3536	c Si c'est la fin, il faut conserver l'etat de redemarrage
3537	c====================================================================
3538	c
3539
3540
3541	IF (lafin) THEN
3542	itau_phy = itau_phy + itap
3543	CALL phyredem ("restartphy.nc")
3544	! open(97,form="unformatted",file="finbin")
3545	! write(97) u_seri,v_seri,t_seri,q_seri
3546	! close(97)
3547	C$OMP MASTER
3548	if (read_climoz >= 1) then
3549	if (is_mpi_root) then
3550	call nf95_close(ncid_climoz)
3551	end if
3552	deallocate(press_climoz) ! pointer
3553	end if
3554	C$OMP END MASTER
3555	ENDIF
3556
3557	! first=.false.
3558
3559	RETURN
3560	END
3561	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3562	IMPLICIT none
3563	c
3564	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3565	c la conservation de l'eau
3566	c
3567	#include "YOMCST.h"
3568	INTEGER klon,klev
3569	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3570	REAL aire(klon)
3571	REAL qtotal, zx, qcheck
3572	INTEGER i, k
3573	c
3574	zx = 0.0
3575	DO i = 1, klon
3576	zx = zx + aire(i)
3577	ENDDO
3578	qtotal = 0.0
3579	DO k = 1, klev
3580	DO i = 1, klon
3581	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3582	. *(paprs(i,k)-paprs(i,k+1))/RG
3583	ENDDO
3584	ENDDO
3585	c
3586	qcheck = qtotal/zx
3587	c
3588	RETURN
3589	END
3590	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3591	IMPLICIT none
3592	c
3593	c Tranformer une variable de la grille physique a
3594	c la grille d'ecriture
3595	c
3596	INTEGER nfield,nlon,iim,jjmp1, jjm
3597	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3598	c
3599	INTEGER i, n, ig
3600	c
3601	jjm = jjmp1 - 1
3602	DO n = 1, nfield
3603	DO i=1,iim
3604	ecrit(i,n) = fi(1,n)
3605	ecrit(i+jjm*iim,n) = fi(nlon,n)
3606	ENDDO
3607	DO ig = 1, nlon - 2
3608	ecrit(iim+ig,n) = fi(1+ig,n)
3609	ENDDO
3610	ENDDO
3611	RETURN
3612	END

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

Download in other formats:

Original Format