Context Navigation

physiq.F @ 1214

Last change on this file since 1214 was 1213, checked in by idelkadi, 15 years ago
Rajout de la possibilite de lire dans *.def les frequences de sorties pour les different fichiers de sortie, sous format : phys_out_filetimesteps= Xmth, Yday, Zhor, Ahor, Bmth (X, Y, Z, A et B entiers) Corrections dans phys_output_mod.F90
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 114.1 KB

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

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

Download in other formats:

Original Format