Context Navigation

physiq.F @ 1201

Last change on this file since 1201 was 1201, checked in by Laurent Fairhead, 16 years ago
Modifications nécessaires a l'inclusion d'un calendrier réaliste. La date courante est calculée dans leapfrog.F et exprimée en Jour Julien (modifié). On en a profité pour faire un peu de ménage dans la gestion des dates du modèle. Dans la physique, on utilise les routines de passages entre calendrier Julien et Gregorien incluses dans IOIPSL pour calculer le nombre de jours écoulés depuis le 1er janvier (pour les conditions aux limites) ou l'equinoxe (pour le calcul de la longitude solaire). Le calcul de l'orbite reprend celui du gcm planétaire (codé par FH) On décide du calendrier à utiliser à l'aide du paramètre calend du run.def. Par défaut celui-ci est à earth_360d LF
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 1201 2009-07-07 14:01:00Z fairhead $
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 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1437	ecrit_hf2mth = ecrit_mth/ecrit_hf
1438
1439	ecrit_hf = ecrit_hf * un_jour
1440	!IM
1441	IF(ecrit_day.LE.1.) THEN
1442	ecrit_day = ecrit_day * un_jour !en secondes
1443	ENDIF
1444	!IM
1445	ecrit_mth = ecrit_mth * un_jour
1446	ecrit_ins = ecrit_ins * un_jour
1447	ecrit_reg = ecrit_reg * un_jour
1448	ecrit_tra = ecrit_tra * un_jour
1449	ecrit_ISCCP = ecrit_ISCCP * un_jour
1450	ecrit_LES = ecrit_LES * un_jour
1451	c
1452	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1453	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1454	. ecrit_hf2mth
1455	cIM 030306 END
1456
1457	capemaxcels = 't_max(X)'
1458	t2mincels = 't_min(X)'
1459	t2maxcels = 't_max(X)'
1460	tinst = 'inst(X)'
1461	tave = 'ave(X)'
1462	cIM cf. AM 081204 BEG
1463	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1464	cIM cf. AM 081204 END
1465	c
1466	c=============================================================
1467	c Initialisation des sorties
1468	c=============================================================
1469
1470	#ifdef CPP_IOIPSL
1471
1472	c$OMP MASTER
1473	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1474	& ctetaSTD,dtime,ok_veget,
1475	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1476	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie)
1477	c$OMP END MASTER
1478	c$OMP BARRIER
1479
1480	#ifdef histISCCP
1481	#include "ini_histISCCP.h"
1482	#endif
1483
1484	#ifdef histmthNMC
1485	#include "ini_histmthNMC.h"
1486	#endif
1487
1488	#include "ini_histday_seri.h"
1489
1490	#include "ini_paramLMDZ_phy.h"
1491
1492	#endif
1493
1494	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1495	date0 = jD_ref
1496	WRITE(,) 'physiq date0 : ',date0
1497	c
1498	c
1499	c
1500	c Prescrire l'ozone dans l'atmosphere
1501	c
1502	c
1503	cc DO i = 1, klon
1504	cc DO k = 1, klev
1505	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1506	cc ENDDO
1507	cc ENDDO
1508	c
1509	IF (config_inca /= 'none') THEN
1510	#ifdef INCA
1511	CALL VTe(VTphysiq)
1512	CALL VTb(VTinca)
1513	! iii = MOD(NINT(xjour),360)
1514	! calday = FLOAT(iii) + jH_cur
1515	calday = FLOAT(days_elapsed) + jH_cur
1516	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1517
1518	CALL chemini(
1519	$ rg,
1520	$ ra,
1521	$ airephy,
1522	$ rlat,
1523	$ rlon,
1524	$ presnivs,
1525	$ calday,
1526	$ klon,
1527	$ nqtot,
1528	$ pdtphys,
1529	$ annee_ref,
1530	$ day_ini)
1531
1532	CALL VTe(VTinca)
1533	CALL VTb(VTphysiq)
1534	#endif
1535	END IF
1536	c
1537	c
1538	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1539	! Nouvelle initialisation pour le rayonnement RRTM
1540	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1541
1542	call iniradia(klon,klev,paprs(1,1:klev+1))
1543
1544	C$omp single
1545	if (read_climoz) then
1546	call open_climoz(ncid_climoz, press_climoz)
1547	END IF
1548	C$omp end single
1549	ENDIF
1550	!
1551	! ************** Fin de IF ( debut ) *************
1552	!
1553	!
1554	! Incrementer le compteur de la physique
1555	!
1556	itap = itap + 1
1557	!
1558	! Update fraction of the sub-surfaces (pctsrf) and
1559	! initialize, where a new fraction has appeared, all variables depending
1560	! on the surface fraction.
1561	!
1562	CALL change_srf_frac(itap, dtime, days_elapsed+1,
1563	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1564
1565	! Tendances bidons pour les processus qui n'affectent pas certaines
1566	! variables.
1567	du0(:,:)=0.
1568	dv0(:,:)=0.
1569	dq0(:,:)=0.
1570	dql0(:,:)=0.
1571	c
1572	c Mettre a zero des variables de sortie (pour securite)
1573	c
1574	DO i = 1, klon
1575	d_ps(i) = 0.0
1576	ENDDO
1577	DO k = 1, klev
1578	DO i = 1, klon
1579	d_t(i,k) = 0.0
1580	d_u(i,k) = 0.0
1581	d_v(i,k) = 0.0
1582	ENDDO
1583	ENDDO
1584	DO iq = 1, nqtot
1585	DO k = 1, klev
1586	DO i = 1, klon
1587	d_qx(i,k,iq) = 0.0
1588	ENDDO
1589	ENDDO
1590	ENDDO
1591	da(:,:)=0.
1592	mp(:,:)=0.
1593	phi(:,:,:)=0.
1594	c
1595	c Ne pas affecter les valeurs entrees de u, v, h, et q
1596	c
1597	DO k = 1, klev
1598	DO i = 1, klon
1599	t_seri(i,k) = t(i,k)
1600	u_seri(i,k) = u(i,k)
1601	v_seri(i,k) = v(i,k)
1602	q_seri(i,k) = qx(i,k,ivap)
1603	ql_seri(i,k) = qx(i,k,iliq)
1604	qs_seri(i,k) = 0.
1605	ENDDO
1606	ENDDO
1607	IF (nqtot.GE.3) THEN
1608	DO iq = 3, nqtot
1609	DO k = 1, klev
1610	DO i = 1, klon
1611	tr_seri(i,k,iq-2) = qx(i,k,iq)
1612	ENDDO
1613	ENDDO
1614	ENDDO
1615	ELSE
1616	DO k = 1, klev
1617	DO i = 1, klon
1618	tr_seri(i,k,1) = 0.0
1619	ENDDO
1620	ENDDO
1621	ENDIF
1622	C
1623	DO i = 1, klon
1624	ztsol(i) = 0.
1625	ENDDO
1626	DO nsrf = 1, nbsrf
1627	DO i = 1, klon
1628	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1629	ENDDO
1630	ENDDO
1631	cIM
1632	IF (ip_ebil_phy.ge.1) THEN
1633	ztit='after dynamic'
1634	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1635	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1636	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1637	C Comme les tendances de la physique sont ajoute dans la dynamique,
1638	C on devrait avoir que la variation d'entalpie par la dynamique
1639	C est egale a la variation de la physique au pas de temps precedent.
1640	C Donc la somme de ces 2 variations devrait etre nulle.
1641	call diagphy(airephy,ztit,ip_ebil_phy
1642	e , zero_v, zero_v, zero_v, zero_v, zero_v
1643	e , zero_v, zero_v, zero_v, ztsol
1644	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1645	s , fs_bound, fq_bound )
1646	END IF
1647
1648	c Diagnostiquer la tendance dynamique
1649	c
1650	IF (ancien_ok) THEN
1651	DO k = 1, klev
1652	DO i = 1, klon
1653	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1654	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1655	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1656	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1657	ENDDO
1658	ENDDO
1659	ELSE
1660	DO k = 1, klev
1661	DO i = 1, klon
1662	d_u_dyn(i,k) = 0.0
1663	d_v_dyn(i,k) = 0.0
1664	d_t_dyn(i,k) = 0.0
1665	d_q_dyn(i,k) = 0.0
1666	ENDDO
1667	ENDDO
1668	ancien_ok = .TRUE.
1669	ENDIF
1670	c
1671	c Ajouter le geopotentiel du sol:
1672	c
1673	DO k = 1, klev
1674	DO i = 1, klon
1675	zphi(i,k) = pphi(i,k) + pphis(i)
1676	ENDDO
1677	ENDDO
1678	c
1679	c Verifier les temperatures
1680	c
1681	cIM BEG
1682	IF (check) THEN
1683	amn=MIN(ftsol(1,is_ter),1000.)
1684	amx=MAX(ftsol(1,is_ter),-1000.)
1685	DO i=2, klon
1686	amn=MIN(ftsol(i,is_ter),amn)
1687	amx=MAX(ftsol(i,is_ter),amx)
1688	ENDDO
1689	c
1690	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1691	ENDIF !(check) THEN
1692	cIM END
1693	c
1694	CALL hgardfou(t_seri,ftsol,'debutphy')
1695	c
1696	cIM BEG
1697	IF (check) THEN
1698	amn=MIN(ftsol(1,is_ter),1000.)
1699	amx=MAX(ftsol(1,is_ter),-1000.)
1700	DO i=2, klon
1701	amn=MIN(ftsol(i,is_ter),amn)
1702	amx=MAX(ftsol(i,is_ter),amx)
1703	ENDDO
1704	c
1705	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1706	ENDIF !(check) THEN
1707	cIM END
1708	c
1709	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1710	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1711	c
1712	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1713	C Once per day, update ozone:
1714	if (read_climoz) then
1715	C Ozone climatology from a NetCDF file
1716	call regr_pr_av(ncid_climoz, "tro3", days_elapsed+1,
1717	& press_climoz,
1718	$ paprs, wo)
1719	! Convert from mole fraction of ozone to column density of ozone in a
1720	! cell, in kDU:
1721	wo = wo * 48. / 29. * zmasse / dobson_u / 1e3
1722	C (By regridding ozone values for LMDZ only once per day, we
1723	C have already neglected the variation of pressure in one
1724	C day. So do not recompute "wo" at each time step even if
1725	C "zmasse" changes a little.)
1726	else
1727	CALL ozonecm(real(days_elapsed+1), rlat, paprs, wo)
1728
1729	end if
1730	ENDIF
1731	c
1732	c Re-evaporer l'eau liquide nuageuse
1733	c
1734	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1735	DO i = 1, klon
1736	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1737	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1738	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1739	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1740	zb = MAX(0.0,ql_seri(i,k))
1741	za = - MAX(0.0,ql_seri(i,k))
1742	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1743	t_seri(i,k) = t_seri(i,k) + za
1744	q_seri(i,k) = q_seri(i,k) + zb
1745	ql_seri(i,k) = 0.0
1746	d_t_eva(i,k) = za
1747	d_q_eva(i,k) = zb
1748	ENDDO
1749	ENDDO
1750	cIM
1751	IF (ip_ebil_phy.ge.2) THEN
1752	ztit='after reevap'
1753	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1754	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1755	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1756	call diagphy(airephy,ztit,ip_ebil_phy
1757	e , zero_v, zero_v, zero_v, zero_v, zero_v
1758	e , zero_v, zero_v, zero_v, ztsol
1759	e , d_h_vcol, d_qt, d_ec
1760	s , fs_bound, fq_bound )
1761	C
1762	END IF
1763
1764	c
1765	c=========================================================================
1766	! Calculs de l'orbite.
1767	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1768	! doit donc etre placé avant radlwsw et pbl_surface
1769
1770	! calcul selon la routine utilisee pour les planetes
1771	if (new_orbit) then
1772	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1773	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1774	! day_since_equinox = (jD_cur) - jD_eq
1775	call solarlong(day_since_equinox, zlongi, dist)
1776	else
1777	! calcul selon la routine utilisee pour l'AR4
1778	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1779	! solarlong0
1780	if (solarlong0<-999.) then
1781	CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
1782	else
1783	zlongi=solarlong0 ! longitude solaire vraie
1784	dist=1. ! distance au soleil / moyenne
1785	endif
1786	endif
1787	if(prt_level.ge.1) &
1788	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1789
1790	! Avec ou sans cycle diurne
1791	IF (cycle_diurne) THEN
1792	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1793	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1794	ELSE
1795	CALL angle(zlongi, rlat, fract, rmu0)
1796	ENDIF
1797
1798	if (mydebug) then
1799	call writefield_phy('u_seri',u_seri,llm)
1800	call writefield_phy('v_seri',v_seri,llm)
1801	call writefield_phy('t_seri',t_seri,llm)
1802	call writefield_phy('q_seri',q_seri,llm)
1803	endif
1804
1805	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1806	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1807	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1808	c turbulent du couche limit.
1809	c
1810	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1811	c ecriture des fihiers hist_XXXX.nc, ces sont :
1812	c qsol, zq2m, s_pblh, s_lcl,
1813	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1814	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1815	c zxrugs, zu10m, zv10m, fder,
1816	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1817	c frugs, agesno, fsollw, fsolsw,
1818	c d_ts, fevap, fluxlat, t2m,
1819	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1820	c
1821	c Certains ne sont pas utiliser du tout :
1822	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1823	c
1824
1825	CALL pbl_surface(
1826	e dtime, date0, itap, days_elapsed+1,
1827	e debut, lafin,
1828	e rlon, rlat, rugoro, rmu0,
1829	e rain_fall, snow_fall, solsw, sollw,
1830	e t_seri, q_seri, u_seri, v_seri,
1831	e pplay, paprs, pctsrf,
1832	+ ftsol, falb1, falb2, u10m, v10m,
1833	s sollwdown, cdragh, cdragm, u1, v1,
1834	s albsol1, albsol2, sens, evap,
1835	s zxtsol, zxfluxlat, zt2m, qsat2m,
1836	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1837	s coefh, slab_wfbils,
1838	d qsol, zq2m, s_pblh, s_lcl,
1839	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1840	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1841	d zxrugs, zu10m, zv10m, fder,
1842	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1843	d frugs, agesno, fsollw, fsolsw,
1844	d d_ts, fevap, fluxlat, t2m,
1845	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1846	- dsens, devap, zxsnow,
1847	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1848
1849
1850	!-----------------------------------------------------------------------------------------
1851	! ajout des tendances de la diffusion turbulente
1852	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1853	!-----------------------------------------------------------------------------------------
1854
1855	if (mydebug) then
1856	call writefield_phy('u_seri',u_seri,llm)
1857	call writefield_phy('v_seri',v_seri,llm)
1858	call writefield_phy('t_seri',t_seri,llm)
1859	call writefield_phy('q_seri',q_seri,llm)
1860	endif
1861
1862
1863	IF (ip_ebil_phy.ge.2) THEN
1864	ztit='after surface_main'
1865	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1866	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1867	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1868	call diagphy(airephy,ztit,ip_ebil_phy
1869	e , zero_v, zero_v, zero_v, zero_v, sens
1870	e , evap , zero_v, zero_v, ztsol
1871	e , d_h_vcol, d_qt, d_ec
1872	s , fs_bound, fq_bound )
1873	END IF
1874
1875	c =================================================================== c
1876	c Calcul de Qsat
1877
1878	DO k = 1, klev
1879	DO i = 1, klon
1880	zx_t = t_seri(i,k)
1881	IF (thermcep) THEN
1882	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1883	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1884	zx_qs = MIN(0.5,zx_qs)
1885	zcor = 1./(1.-retv*zx_qs)
1886	zx_qs = zx_qs*zcor
1887	ELSE
1888	IF (zx_t.LT.t_coup) THEN
1889	zx_qs = qsats(zx_t)/pplay(i,k)
1890	ELSE
1891	zx_qs = qsatl(zx_t)/pplay(i,k)
1892	ENDIF
1893	ENDIF
1894	zqsat(i,k)=zx_qs
1895	ENDDO
1896	ENDDO
1897
1898	if (prt_level.ge.1) then
1899	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1900	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1901	endif
1902	c
1903	c Appeler la convection (au choix)
1904	c
1905	DO k = 1, klev
1906	DO i = 1, klon
1907	conv_q(i,k) = d_q_dyn(i,k)
1908	. + d_q_vdf(i,k)/dtime
1909	conv_t(i,k) = d_t_dyn(i,k)
1910	. + d_t_vdf(i,k)/dtime
1911	ENDDO
1912	ENDDO
1913	IF (check) THEN
1914	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1915	WRITE(lunout,*) "avantcon=", za
1916	ENDIF
1917	zx_ajustq = .FALSE.
1918	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1919	IF (zx_ajustq) THEN
1920	DO i = 1, klon
1921	z_avant(i) = 0.0
1922	ENDDO
1923	DO k = 1, klev
1924	DO i = 1, klon
1925	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1926	. *(paprs(i,k)-paprs(i,k+1))/RG
1927	ENDDO
1928	ENDDO
1929	ENDIF
1930
1931	c Calcule de vitesse verticale a partir de flux de masse verticale
1932	DO k = 1, klev
1933	DO i = 1, klon
1934	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1935	END DO
1936	END DO
1937
1938	IF (iflag_con.EQ.1) THEN
1939	stop'reactiver le call conlmd dans physiq.F'
1940	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1941	c . d_t_con, d_q_con,
1942	c . rain_con, snow_con, ibas_con, itop_con)
1943	ELSE IF (iflag_con.EQ.2) THEN
1944	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1945	e conv_t, conv_q, -evap, omega,
1946	s d_t_con, d_q_con, rain_con, snow_con,
1947	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1948	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1949	d_u_con = 0.
1950	d_v_con = 0.
1951
1952	WHERE (rain_con < 0.) rain_con = 0.
1953	WHERE (snow_con < 0.) snow_con = 0.
1954	DO i = 1, klon
1955	ibas_con(i) = klev+1 - kcbot(i)
1956	itop_con(i) = klev+1 - kctop(i)
1957	ENDDO
1958	ELSE IF (iflag_con.GE.3) THEN
1959	c nb of tracers for the KE convection:
1960	c MAF la partie traceurs est faite dans phytrac
1961	c on met ntra=1 pour limiter les appels mais on peut
1962	c supprimer les calculs / ftra.
1963	ntra = 1
1964
1965	c=====================================================================================
1966	cajout pour la parametrisation des poches froides:
1967	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1968	do k=1,klev
1969	do i=1,klon
1970	if (iflag_wake.eq.1) then
1971	t_wake(i,k) = t_seri(i,k)
1972	. +(1-wake_s(i))*wake_deltat(i,k)
1973	q_wake(i,k) = q_seri(i,k)
1974	. +(1-wake_s(i))*wake_deltaq(i,k)
1975	t_undi(i,k) = t_seri(i,k)
1976	. -wake_s(i)*wake_deltat(i,k)
1977	q_undi(i,k) = q_seri(i,k)
1978	. -wake_s(i)*wake_deltaq(i,k)
1979	else
1980	t_wake(i,k) = t_seri(i,k)
1981	q_wake(i,k) = q_seri(i,k)
1982	t_undi(i,k) = t_seri(i,k)
1983	q_undi(i,k) = q_seri(i,k)
1984	endif
1985	enddo
1986	enddo
1987
1988	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
1989	cc-- pour le soulevement des particules dans le modele convectif
1990	c
1991	do i = 1,klon
1992	ALE(i) = 0.
1993	ALP(i) = 0.
1994	enddo
1995	c
1996	ccalcul de ale_wake et alp_wake
1997	do i = 1,klon
1998	if (iflag_wake.eq.1) then
1999	ale_wake(i) = 0.5wake_cstar(i)*2
2000	alp_wake(i) = wake_fip(i)
2001	else
2002	ale_wake(i) = 0.
2003	alp_wake(i) = 0.
2004	endif
2005	enddo
2006	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2007	cdans le thermique sinon
2008	if (iflag_coupl.eq.0) then
2009	if (debut) print*,'ALE et ALP imposes'
2010	do i = 1,klon
2011	con ne couple que ale
2012	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2013	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2014	con ne couple que alp
2015	c ALP(i) = alp_wake(i) + Alp_bl(i)
2016	ALP(i) = alp_wake(i) + alp_bl_prescr
2017	enddo
2018	else
2019	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2020	do i = 1,klon
2021	ALE(i) = max(ale_wake(i),Ale_bl(i))
2022	ALP(i) = alp_wake(i) + Alp_bl(i)
2023	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2024	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2025	enddo
2026	endif
2027	do i=1,klon
2028	if (alp(i)>alp_max) then
2029	IF(prt_level>9)WRITE(lunout,*) &
2030	& 'WARNING SUPER ALP (seuil=',alp_max,
2031	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2032	alp(i)=alp_max
2033	endif
2034	if (ale(i)>ale_max) then
2035	IF(prt_level>9)WRITE(lunout,*) &
2036	& 'WARNING SUPER ALE (seuil=',ale_max,
2037	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2038	ale(i)=ale_max
2039	endif
2040	enddo
2041
2042	cfin calcul ale et alp
2043	c=================================================================================================
2044
2045
2046	c sb, oct02:
2047	c Schema de convection modularise et vectorise:
2048	c (driver commun aux versions 3 et 4)
2049	c
2050	IF (ok_cvl) THEN ! new driver for convectL
2051
2052	CALL concvl (iflag_con,iflag_clos,
2053	. dtime,paprs,pplay,t_undi,q_undi,
2054	. t_wake,q_wake,wake_s,
2055	. u_seri,v_seri,tr_seri,nbtr,
2056	. ALE,ALP,
2057	. ema_work1,ema_work2,
2058	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2059	. rain_con, snow_con, ibas_con, itop_con, sigd,
2060	. upwd,dnwd,dnwd0,
2061	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2062	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2063	. pmflxr,pmflxs,da,phi,mp,
2064	. ftd,fqd,lalim_conv,wght_th)
2065
2066	cIM begin
2067	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2068	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2069	cIM end
2070	cIM cf. FH
2071	clwcon0=qcondc
2072	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2073
2074	ELSE ! ok_cvl
2075	c MAF conema3 ne contient pas les traceurs
2076	CALL conema3 (dtime,
2077	. paprs,pplay,t_seri,q_seri,
2078	. u_seri,v_seri,tr_seri,ntra,
2079	. ema_work1,ema_work2,
2080	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2081	. rain_con, snow_con, ibas_con, itop_con,
2082	. upwd,dnwd,dnwd0,bas,top,
2083	. Ma,cape,tvp,rflag,
2084	. pbase
2085	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2086	. ,clwcon0)
2087
2088	ENDIF ! ok_cvl
2089
2090	c
2091	c Correction precip
2092	rain_con = rain_con * cvl_corr
2093	snow_con = snow_con * cvl_corr
2094	c
2095
2096	IF (.NOT. ok_gust) THEN
2097	do i = 1, klon
2098	wd(i)=0.0
2099	enddo
2100	ENDIF
2101
2102	c =================================================================== c
2103	c Calcul des proprietes des nuages convectifs
2104	c
2105
2106	c calcul des proprietes des nuages convectifs
2107	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2108	call clouds_gno
2109	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2110
2111	c =================================================================== c
2112
2113	DO i = 1, klon
2114	ema_pcb(i) = pbase(i)
2115	ENDDO
2116	DO i = 1, klon
2117
2118	! L'idicage de itop_con peut cacher un pb potentiel
2119	! FH sous la dictee de JYG, CR
2120	ema_pct(i) = paprs(i,itop_con(i)+1)
2121
2122	if (itop_con(i).gt.klev-3) then
2123	if(prt_level >= 9) then
2124	write(lunout,*)'La convection monte trop haut '
2125	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2126	endif
2127	endif
2128	ENDDO
2129	DO i = 1, klon
2130	ema_cbmf(i) = ema_workcbmf(i)
2131	ENDDO
2132	ELSE IF (iflag_con.eq.0) THEN
2133	write(lunout,*) 'On n appelle pas la convection'
2134	clwcon0=0.
2135	rnebcon0=0.
2136	d_t_con=0.
2137	d_q_con=0.
2138	d_u_con=0.
2139	d_v_con=0.
2140	rain_con=0.
2141	snow_con=0.
2142	bas=1
2143	top=1
2144	ELSE
2145	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2146	CALL abort
2147	ENDIF
2148
2149	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2150	c . d_u_con, d_v_con)
2151
2152	!-----------------------------------------------------------------------------------------
2153	! ajout des tendances de la diffusion turbulente
2154	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2155	!-----------------------------------------------------------------------------------------
2156
2157	if (mydebug) then
2158	call writefield_phy('u_seri',u_seri,llm)
2159	call writefield_phy('v_seri',v_seri,llm)
2160	call writefield_phy('t_seri',t_seri,llm)
2161	call writefield_phy('q_seri',q_seri,llm)
2162	endif
2163
2164	cIM
2165	IF (ip_ebil_phy.ge.2) THEN
2166	ztit='after convect'
2167	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2168	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2169	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2170	call diagphy(airephy,ztit,ip_ebil_phy
2171	e , zero_v, zero_v, zero_v, zero_v, zero_v
2172	e , zero_v, rain_con, snow_con, ztsol
2173	e , d_h_vcol, d_qt, d_ec
2174	s , fs_bound, fq_bound )
2175	END IF
2176	C
2177	IF (check) THEN
2178	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2179	WRITE(lunout,*)"aprescon=", za
2180	zx_t = 0.0
2181	za = 0.0
2182	DO i = 1, klon
2183	za = za + airephy(i)/FLOAT(klon)
2184	zx_t = zx_t + (rain_con(i)+
2185	. snow_con(i))*airephy(i)/FLOAT(klon)
2186	ENDDO
2187	zx_t = zx_t/za*dtime
2188	WRITE(lunout,*)"Precip=", zx_t
2189	ENDIF
2190	IF (zx_ajustq) THEN
2191	DO i = 1, klon
2192	z_apres(i) = 0.0
2193	ENDDO
2194	DO k = 1, klev
2195	DO i = 1, klon
2196	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2197	. *(paprs(i,k)-paprs(i,k+1))/RG
2198	ENDDO
2199	ENDDO
2200	DO i = 1, klon
2201	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2202	. /z_apres(i)
2203	ENDDO
2204	DO k = 1, klev
2205	DO i = 1, klon
2206	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2207	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2208	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2209	ENDIF
2210	ENDDO
2211	ENDDO
2212	ENDIF
2213	zx_ajustq=.FALSE.
2214
2215	c
2216	c=============================================================================
2217	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2218	cpour la couche limite diffuse pour l instant
2219	c
2220	if (iflag_wake.eq.1) then
2221	DO k=1,klev
2222	DO i=1,klon
2223	dt_dwn(i,k) = ftd(i,k)
2224	wdt_PBL(i,k) = 0.
2225	dq_dwn(i,k) = fqd(i,k)
2226	wdq_PBL(i,k) = 0.
2227	M_dwn(i,k) = dnwd0(i,k)
2228	M_up(i,k) = upwd(i,k)
2229	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2230	udt_PBL(i,k) = 0.
2231	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2232	udq_PBL(i,k) = 0.
2233	ENDDO
2234	ENDDO
2235	c
2236	ccalcul caracteristiques de la poche froide
2237	call calWAKE (paprs,pplay,dtime
2238	: ,t_seri,q_seri,omega
2239	: ,dt_dwn,dq_dwn,M_dwn,M_up
2240	: ,dt_a,dq_a,sigd
2241	: ,wdt_PBL,wdq_PBL
2242	: ,udt_PBL,udq_PBL
2243	o ,wake_deltat,wake_deltaq,wake_dth
2244	o ,wake_h,wake_s,wake_dens
2245	o ,wake_pe,wake_fip,wake_gfl
2246	o ,dt_wake,dq_wake
2247	o ,wake_k, t_undi,q_undi
2248	o ,wake_omgbdth,wake_dp_omgb
2249	o ,wake_dtKE,wake_dqKE
2250	o ,wake_dtPBL,wake_dqPBL
2251	o ,wake_omg,wake_dp_deltomg
2252	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2253	o ,wake_ddeltat,wake_ddeltaq)
2254	c
2255	!-----------------------------------------------------------------------------------------
2256	! ajout des tendances des poches froides
2257	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2258	! coherent avec les autres d_t_...
2259	d_t_wake(:,:)=dt_wake(:,:)*dtime
2260	d_q_wake(:,:)=dq_wake(:,:)*dtime
2261	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2262	!-----------------------------------------------------------------------------------------
2263
2264	endif
2265	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2266	c
2267	c===================================================================
2268	c Convection seche (thermiques ou ajustement)
2269	c===================================================================
2270	c
2271	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2272	s ,seuil_inversion,weak_inversion,dthmin)
2273
2274
2275
2276	d_t_ajsb(:,:)=0.
2277	d_q_ajsb(:,:)=0.
2278	d_t_ajs(:,:)=0.
2279	d_u_ajs(:,:)=0.
2280	d_v_ajs(:,:)=0.
2281	d_q_ajs(:,:)=0.
2282	clwcon0th(:,:)=0.
2283	c
2284	fm_therm(:,:)=0.
2285	entr_therm(:,:)=0.
2286	detr_therm(:,:)=0.
2287	c
2288	IF(prt_level>9)WRITE(lunout,*)
2289	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2290	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2291	if(iflag_thermals.lt.0) then
2292	c Rien
2293	c ====
2294	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2295
2296
2297	else
2298
2299	c Thermiques
2300	c ==========
2301	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2302	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2303
2304
2305	if (iflag_thermals.gt.1) then
2306	call calltherm(pdtphys
2307	s ,pplay,paprs,pphi,weak_inversion
2308	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2309	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2310	s ,fm_therm,entr_therm,detr_therm
2311	s ,zqasc,clwcon0th,lmax_th,ratqscth
2312	s ,ratqsdiff,zqsatth
2313	con rajoute ale et alp, et les caracteristiques de la couche alim
2314	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
2315	endif
2316
2317
2318	c Ajustement sec
2319	c ==============
2320
2321	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2322	! a partir du sommet des thermiques.
2323	! Dans le cas contraire, on demarre au niveau 1.
2324
2325	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2326
2327	if(iflag_thermals.eq.0) then
2328	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2329	limbas(:)=1
2330	else
2331	limbas(:)=lmax_th(:)
2332	endif
2333
2334	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2335	! pour des test de convergence numerique.
2336	! Le nouveau ajsec est a priori mieux, meme pour le cas
2337	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2338	! non nulles numeriquement pour des mailles non concernees.
2339
2340	if (iflag_thermals.eq.0) then
2341	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2342	s , d_t_ajsb, d_q_ajsb)
2343	else
2344	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2345	s , d_t_ajsb, d_q_ajsb)
2346	endif
2347
2348	!-----------------------------------------------------------------------------------------
2349	! ajout des tendances de l'ajustement sec ou des thermiques
2350	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2351	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2352	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2353
2354	!-----------------------------------------------------------------------------------------
2355
2356	endif
2357
2358	endif
2359	c
2360	c===================================================================
2361	cIM
2362	IF (ip_ebil_phy.ge.2) THEN
2363	ztit='after dry_adjust'
2364	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2365	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2366	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2367	END IF
2368
2369
2370	c-------------------------------------------------------------------------
2371	c Caclul des ratqs
2372	c-------------------------------------------------------------------------
2373
2374	c print*,'calcul des ratqs'
2375	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2376	c ----------------
2377	c on ecrase le tableau ratqsc calcule par clouds_gno
2378	if (iflag_cldcon.eq.1) then
2379	do k=1,klev
2380	do i=1,klon
2381	if(ptconv(i,k)) then
2382	ratqsc(i,k)=ratqsbas
2383	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2384	else
2385	ratqsc(i,k)=0.
2386	endif
2387	enddo
2388	enddo
2389
2390	c-----------------------------------------------------------------------
2391	c par nversion de la fonction log normale
2392	c-----------------------------------------------------------------------
2393	else if (iflag_cldcon.eq.4) then
2394	ptconvth(:,:)=.false.
2395	ratqsc(:,:)=0.
2396	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2397	call clouds_gno
2398	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2399	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2400
2401	c-----------------------------------------------------------------------
2402	c par calcul direct de l'ecart-type
2403	c-----------------------------------------------------------------------
2404
2405	else if (iflag_cldcon>=5) then
2406	wmax_th(:)=0.
2407	zmax_th(:)=0.
2408	do k=1,klev
2409	do i=1,klon
2410	wmax_th(i)=max(wmax_th(i),zw2(i,k))
2411	if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
2412	enddo
2413	enddo
2414	tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
2415	print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
2416
2417	c On impose que l'air autour de la fraction couverte par le thermique
2418	c plus son air detraine durant tau_overturning_th soit superieur
2419	c a 0.1 q_seri
2420	zz=0.1
2421	do k=1,klev
2422	do i=1,klon
2423	lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
2424	s tau_overturning_th(i)*detr_therm(i,k)
2425	s *rg/(paprs(i,k)-paprs(i,k+1))
2426	znum=(1.-zz)*q_seri(i,k)
2427	zden=zqasc(i,k)-zz*q_seri(i,k)
2428	if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
2429	lambda_th(i,k)=min(lambda_th(i,k),0.9)
2430	enddo
2431	enddo
2432
2433	if(iflag_cldcon==5) then
2434	do k=1,klev
2435	do i=1,klon
2436	ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
2437	s abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
2438	enddo
2439	enddo
2440	else if(iflag_cldcon==6) then
2441	do k=1,klev
2442	do i=1,klon
2443	ratqsc(i,k)=sqrt(lambda_th(i,k))*
2444	s (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
2445	enddo
2446	enddo
2447	endif
2448
2449	endif
2450
2451	c ratqs stables
2452	c -------------
2453
2454	if (iflag_ratqs.eq.0) then
2455
2456	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2457	do k=1,klev
2458	do i=1, klon
2459	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2460	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2461	enddo
2462	enddo
2463
2464	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2465	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2466	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2467	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2468	! Il s'agit de differents tests dans la phase de reglage du modele
2469	! avec thermiques.
2470
2471	else if (iflag_ratqs.eq.1) then
2472
2473	do k=1,klev
2474	do i=1, klon
2475	if (pplay(i,k).ge.60000.) then
2476	ratqss(i,k)=ratqsbas
2477	else if ((pplay(i,k).ge.30000.).and.
2478	s (pplay(i,k).lt.60000.)) then
2479	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2480	s (60000.-pplay(i,k))/(60000.-30000.)
2481	else
2482	ratqss(i,k)=ratqshaut
2483	endif
2484	enddo
2485	enddo
2486
2487	else
2488
2489	do k=1,klev
2490	do i=1, klon
2491	if (pplay(i,k).ge.60000.) then
2492	ratqss(i,k)=ratqsbas
2493	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2494	else if ((pplay(i,k).ge.30000.).and.
2495	s (pplay(i,k).lt.60000.)) then
2496	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2497	s (60000.-pplay(i,k))/(60000.-30000.)
2498	else
2499	ratqss(i,k)=ratqshaut
2500	endif
2501	enddo
2502	enddo
2503	endif
2504
2505
2506
2507
2508	c ratqs final
2509	c -----------
2510
2511	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2512	s .or.iflag_cldcon.ge.4) then
2513
2514	! On ajoute une constante au ratqsc*2 pour tenir compte de
2515	! fluctuations turbulentes de petite echelle
2516
2517	do k=1,klev
2518	do i=1,klon
2519	if ((fm_therm(i,k).gt.1.e-10)) then
2520	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2521	endif
2522	enddo
2523	enddo
2524
2525	! les ratqs sont une conbinaison de ratqss et ratqsc
2526	! 1800s, en dur pour le moment, est le temps de
2527	! relaxation des ratqs
2528
2529	facteur=exp(-pdtphys/1800.)
2530
2531	print*,'WARNING ratqs a revoir '
2532	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2533	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2534
2535	else
2536	! on ne prend que le ratqs stable pour fisrtilp
2537	ratqs(:,:)=ratqss(:,:)
2538	endif
2539
2540
2541	c
2542	c Appeler le processus de condensation a grande echelle
2543	c et le processus de precipitation
2544	c-------------------------------------------------------------------------
2545	CALL fisrtilp(dtime,paprs,pplay,
2546	. t_seri, q_seri,ptconv,ratqs,
2547	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2548	. rain_lsc, snow_lsc,
2549	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2550	. frac_impa, frac_nucl,
2551	. prfl, psfl, rhcl)
2552
2553	WHERE (rain_lsc < 0) rain_lsc = 0.
2554	WHERE (snow_lsc < 0) snow_lsc = 0.
2555	!-----------------------------------------------------------------------------------------
2556	! ajout des tendances de la diffusion turbulente
2557	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2558	!-----------------------------------------------------------------------------------------
2559	DO k = 1, klev
2560	DO i = 1, klon
2561	cldfra(i,k) = rneb(i,k)
2562	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2563	ENDDO
2564	ENDDO
2565	IF (check) THEN
2566	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2567	WRITE(lunout,*)"apresilp=", za
2568	zx_t = 0.0
2569	za = 0.0
2570	DO i = 1, klon
2571	za = za + airephy(i)/FLOAT(klon)
2572	zx_t = zx_t + (rain_lsc(i)
2573	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2574	ENDDO
2575	zx_t = zx_t/za*dtime
2576	WRITE(lunout,*)"Precip=", zx_t
2577	ENDIF
2578	cIM
2579	IF (ip_ebil_phy.ge.2) THEN
2580	ztit='after fisrt'
2581	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2582	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2583	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2584	call diagphy(airephy,ztit,ip_ebil_phy
2585	e , zero_v, zero_v, zero_v, zero_v, zero_v
2586	e , zero_v, rain_lsc, snow_lsc, ztsol
2587	e , d_h_vcol, d_qt, d_ec
2588	s , fs_bound, fq_bound )
2589	END IF
2590
2591	if (mydebug) then
2592	call writefield_phy('u_seri',u_seri,llm)
2593	call writefield_phy('v_seri',v_seri,llm)
2594	call writefield_phy('t_seri',t_seri,llm)
2595	call writefield_phy('q_seri',q_seri,llm)
2596	endif
2597
2598	c
2599	c-------------------------------------------------------------------
2600	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2601	c-------------------------------------------------------------------
2602
2603	c 1. NUAGES CONVECTIFS
2604	c
2605	cIM cf FH
2606	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2607	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2608	snow_tiedtke=0.
2609	c print*,'avant calcul de la pseudo precip '
2610	c print*,'iflag_cldcon',iflag_cldcon
2611	if (iflag_cldcon.eq.-1) then
2612	rain_tiedtke=rain_con
2613	else
2614	c print*,'calcul de la pseudo precip '
2615	rain_tiedtke=0.
2616	c print*,'calcul de la pseudo precip 0'
2617	do k=1,klev
2618	do i=1,klon
2619	if (d_q_con(i,k).lt.0.) then
2620	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2621	s *(paprs(i,k)-paprs(i,k+1))/rg
2622	endif
2623	enddo
2624	enddo
2625	endif
2626	c
2627	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2628	c
2629
2630	c Nuages diagnostiques pour Tiedtke
2631	CALL diagcld1(paprs,pplay,
2632	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2633	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2634	. diafra,dialiq)
2635	DO k = 1, klev
2636	DO i = 1, klon
2637	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2638	cldliq(i,k) = dialiq(i,k)
2639	cldfra(i,k) = diafra(i,k)
2640	ENDIF
2641	ENDDO
2642	ENDDO
2643
2644	ELSE IF (iflag_cldcon.ge.3) THEN
2645	c On prend pour les nuages convectifs le max du calcul de la
2646	c convection et du calcul du pas de temps precedent diminue d'un facteur
2647	c facttemps
2648	facteur = pdtphys *facttemps
2649	do k=1,klev
2650	do i=1,klon
2651	rnebcon(i,k)=rnebcon(i,k)*facteur
2652	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2653	s then
2654	rnebcon(i,k)=rnebcon0(i,k)
2655	clwcon(i,k)=clwcon0(i,k)
2656	endif
2657	enddo
2658	enddo
2659
2660	c
2661	cjq - introduce the aerosol direct and first indirect radiative forcings
2662	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2663	IF (ok_ade.OR.ok_aie) THEN
2664	IF (.NOT. aerosol_couple)
2665	& CALL readaerosol_optic(
2666	& debut, new_aod, flag_aerosol, jD_cur-jD_ref, pdtphys,
2667	& pplay, paprs, t_seri, rhcl,
2668	& mass_solu_aero, mass_solu_aero_pi,
2669	& tau_aero, piz_aero, cg_aero,
2670	& tausum_aero, tau3d_aero)
2671	ELSE
2672	tau_aero(:,:,:,:) = 0.
2673	piz_aero(:,:,:,:) = 0.
2674	cg_aero(:,:,:,:) = 0.
2675	ENDIF
2676
2677	cIM calcul nuages par le simulateur ISCCP
2678	c
2679	#ifdef histISCCP
2680	IF (ok_isccp) THEN
2681	c
2682	cIM lecture invtau, tautab des fichiers formattes
2683	c
2684	IF (debut) THEN
2685	c$OMP MASTER
2686	c
2687	open(99,file='tautab.formatted', FORM='FORMATTED')
2688	read(99,'(f30.20)') tautab_omp
2689	close(99)
2690	c
2691	open(99,file='invtau.formatted',form='FORMATTED')
2692	read(99,'(i10)') invtau_omp
2693
2694	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
2695	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2696
2697	close(99)
2698	c$OMP END MASTER
2699	c$OMP BARRIER
2700	tautab=tautab_omp
2701	invtau=invtau_omp
2702	c
2703	ENDIF !debut
2704	c
2705	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2706	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2707	#include "calcul_simulISCCP.h"
2708	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2709	ENDIF !ok_isccp
2710	#endif
2711
2712	c On prend la somme des fractions nuageuses et des contenus en eau
2713	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2714	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2715
2716	ENDIF
2717	c
2718	c 2. NUAGES STARTIFORMES
2719	c
2720	IF (ok_stratus) THEN
2721	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2722	DO k = 1, klev
2723	DO i = 1, klon
2724	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2725	cldliq(i,k) = dialiq(i,k)
2726	cldfra(i,k) = diafra(i,k)
2727	ENDIF
2728	ENDDO
2729	ENDDO
2730	ENDIF
2731	c
2732	c Precipitation totale
2733	c
2734	DO i = 1, klon
2735	rain_fall(i) = rain_con(i) + rain_lsc(i)
2736	snow_fall(i) = snow_con(i) + snow_lsc(i)
2737	ENDDO
2738	cIM
2739	IF (ip_ebil_phy.ge.2) THEN
2740	ztit="after diagcld"
2741	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2742	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2743	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2744	END IF
2745	c
2746	c Calculer l'humidite relative pour diagnostique
2747	c
2748	DO k = 1, klev
2749	DO i = 1, klon
2750	zx_t = t_seri(i,k)
2751	IF (thermcep) THEN
2752	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2753	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2754	zx_qs = MIN(0.5,zx_qs)
2755	zcor = 1./(1.-retv*zx_qs)
2756	zx_qs = zx_qs*zcor
2757	ELSE
2758	IF (zx_t.LT.t_coup) THEN
2759	zx_qs = qsats(zx_t)/pplay(i,k)
2760	ELSE
2761	zx_qs = qsatl(zx_t)/pplay(i,k)
2762	ENDIF
2763	ENDIF
2764	zx_rh(i,k) = q_seri(i,k)/zx_qs
2765	zqsat(i,k)=zx_qs
2766	ENDDO
2767	ENDDO
2768
2769	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2770	c equivalente a 2m (tpote) pour diagnostique
2771	c
2772	DO i = 1, klon
2773	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2774	IF (thermcep) THEN
2775	IF(zt2m(i).LT.RTT) then
2776	Lheat=RLSTT
2777	ELSE
2778	Lheat=RLVTT
2779	ENDIF
2780	ELSE
2781	IF (zt2m(i).LT.RTT) THEN
2782	Lheat=RLSTT
2783	ELSE
2784	Lheat=RLVTT
2785	ENDIF
2786	ENDIF
2787	tpote(i) = tpot(i)*
2788	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2789	ENDDO
2790
2791	IF (config_inca /= 'none') THEN
2792	#ifdef INCA
2793	CALL VTe(VTphysiq)
2794	CALL VTb(VTinca)
2795	calday = FLOAT(days_elapsed + 1) + jH_cur
2796
2797	IF (config_inca == 'aero') THEN
2798	CALL AEROSOL_METEO_CALC(
2799	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2800	$ prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
2801	END IF
2802
2803	zxsnow_dummy(:) = 0.0
2804
2805	CALL chemhook_begin (calday,
2806	$ days_elapsed,
2807	$ jH_cur,
2808	$ pctsrf(1,1),
2809	$ rlat,
2810	$ rlon,
2811	$ airephy,
2812	$ paprs,
2813	$ pplay,
2814	$ coefh,
2815	$ pphi,
2816	$ t_seri,
2817	$ u,
2818	$ v,
2819	$ wo,
2820	$ q_seri,
2821	$ zxtsol,
2822	$ zxsnow_dummy,
2823	$ solsw,
2824	$ albsol1,
2825	$ rain_fall,
2826	$ snow_fall,
2827	$ itop_con,
2828	$ ibas_con,
2829	$ cldfra,
2830	$ iim,
2831	$ jjm,
2832	$ tr_seri,
2833	$ ftsol,
2834	$ paprs,
2835	$ cdragh,
2836	$ cdragm,
2837	$ pctsrf,
2838	$ pdtphys,
2839	$ itap)
2840
2841	CALL VTe(VTinca)
2842	CALL VTb(VTphysiq)
2843	#endif
2844	END IF !config_inca /= 'none'
2845	c
2846	c Calculer les parametres optiques des nuages et quelques
2847	c parametres pour diagnostiques:
2848	c
2849
2850	IF (aerosol_couple) THEN
2851	mass_solu_aero(:,:) = ccm(:,:,1)
2852	mass_solu_aero_pi(:,:) = ccm(:,:,2)
2853	END IF
2854
2855	if (ok_newmicro) then
2856	CALL newmicro (paprs, pplay,ok_newmicro,
2857	. t_seri, cldliq, cldfra, cldtau, cldemi,
2858	. cldh, cldl, cldm, cldt, cldq,
2859	. flwp, fiwp, flwc, fiwc,
2860	e ok_aie,
2861	e mass_solu_aero, mass_solu_aero_pi,
2862	e bl95_b0, bl95_b1,
2863	s cldtaupi, re, fl)
2864	else
2865	CALL nuage (paprs, pplay,
2866	. t_seri, cldliq, cldfra, cldtau, cldemi,
2867	. cldh, cldl, cldm, cldt, cldq,
2868	e ok_aie,
2869	e mass_solu_aero, mass_solu_aero_pi,
2870	e bl95_b0, bl95_b1,
2871	s cldtaupi, re, fl)
2872
2873	endif
2874	c
2875	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2876	c
2877	IF (MOD(itaprad,radpas).EQ.0) THEN
2878
2879	DO i = 1, klon
2880	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2881	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2882	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2883	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2884	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2885	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2886	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2887	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2888	ENDDO
2889
2890	if (mydebug) then
2891	call writefield_phy('u_seri',u_seri,llm)
2892	call writefield_phy('v_seri',v_seri,llm)
2893	call writefield_phy('t_seri',t_seri,llm)
2894	call writefield_phy('q_seri',q_seri,llm)
2895	endif
2896
2897	IF (aerosol_couple) THEN
2898	#ifdef INCA
2899	CALL radlwsw_inca
2900	e (kdlon,kflev,dist, rmu0, fract, solaire,
2901	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2902	e wo,
2903	e cldfra, cldemi, cldtau,
2904	s heat,heat0,cool,cool0,radsol,albpla,
2905	s topsw,toplw,solsw,sollw,
2906	s sollwdown,
2907	s topsw0,toplw0,solsw0,sollw0,
2908	s lwdn0, lwdn, lwup0, lwup,
2909	s swdn0, swdn, swup0, swup,
2910	e ok_ade, ok_aie,
2911	e tau_aero, piz_aero, cg_aero,
2912	s topswad_aero, solswad_aero,
2913	s topswad0_aero, solswad0_aero,
2914	s topsw_aero, topsw0_aero,
2915	s solsw_aero, solsw0_aero,
2916	e cldtaupi,
2917	s topswai_aero, solswai_aero)
2918
2919	#endif
2920	ELSE
2921
2922	CALL radlwsw
2923	e (dist, rmu0, fract,
2924	e paprs, pplay,zxtsol,albsol1, albsol2,
2925	e t_seri,q_seri,wo,
2926	e cldfra, cldemi, cldtau,
2927	e ok_ade, ok_aie,
2928	e tau_aero, piz_aero, cg_aero,
2929	e cldtaupi,new_aod,
2930	e zqsat, flwc, fiwc,
2931	s heat,heat0,cool,cool0,radsol,albpla,
2932	s topsw,toplw,solsw,sollw,
2933	s sollwdown,
2934	s topsw0,toplw0,solsw0,sollw0,
2935	s lwdn0, lwdn, lwup0, lwup,
2936	s swdn0, swdn, swup0, swup,
2937	s topswad_aero, solswad_aero,
2938	s topswai_aero, solswai_aero,
2939	o topswad0_aero, solswad0_aero,
2940	o topsw_aero, topsw0_aero,
2941	o solsw_aero, solsw0_aero)
2942
2943
2944	ENDIF ! aerosol_couple
2945	itaprad = 0
2946	ENDIF ! MOD(itaprad,radpas)
2947	itaprad = itaprad + 1
2948
2949	if (iflag_radia.eq.0) then
2950	print *,'--------------------------------------------------'
2951	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
2952	print *,'>>>> heat et cool mis a zero '
2953	print *,'--------------------------------------------------'
2954	heat=0.
2955	cool=0.
2956	endif
2957
2958	c
2959	c Ajouter la tendance des rayonnements (tous les pas)
2960	c
2961	DO k = 1, klev
2962	DO i = 1, klon
2963	t_seri(i,k) = t_seri(i,k)
2964	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
2965	ENDDO
2966	ENDDO
2967	c
2968	if (mydebug) then
2969	call writefield_phy('u_seri',u_seri,llm)
2970	call writefield_phy('v_seri',v_seri,llm)
2971	call writefield_phy('t_seri',t_seri,llm)
2972	call writefield_phy('q_seri',q_seri,llm)
2973	endif
2974
2975	cIM
2976	IF (ip_ebil_phy.ge.2) THEN
2977	ztit='after rad'
2978	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2979	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2980	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2981	call diagphy(airephy,ztit,ip_ebil_phy
2982	e , topsw, toplw, solsw, sollw, zero_v
2983	e , zero_v, zero_v, zero_v, ztsol
2984	e , d_h_vcol, d_qt, d_ec
2985	s , fs_bound, fq_bound )
2986	END IF
2987	c
2988	c
2989	c Calculer l'hydrologie de la surface
2990	c
2991	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2992	c . agesno, ftsol,fqsurf,fsnow, ruis)
2993	c
2994
2995	c
2996	c Calculer le bilan du sol et la derive de temperature (couplage)
2997	c
2998	DO i = 1, klon
2999	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3000	c a la demande de JLD
3001	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3002	ENDDO
3003	c
3004	cmoddeblott(jan95)
3005	c Appeler le programme de parametrisation de l'orographie
3006	c a l'echelle sous-maille:
3007	c
3008	IF (ok_orodr) THEN
3009	c
3010	c selection des points pour lesquels le shema est actif:
3011	igwd=0
3012	DO i=1,klon
3013	itest(i)=0
3014	c IF ((zstd(i).gt.10.0)) THEN
3015	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3016	itest(i)=1
3017	igwd=igwd+1
3018	idx(igwd)=i
3019	ENDIF
3020	ENDDO
3021	c igwdim=MAX(1,igwd)
3022	c
3023	IF (ok_strato) THEN
3024
3025	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3026	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3027	e igwd,idx,itest,
3028	e t_seri, u_seri, v_seri,
3029	s zulow, zvlow, zustrdr, zvstrdr,
3030	s d_t_oro, d_u_oro, d_v_oro)
3031
3032	ELSE
3033	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3034	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3035	e igwd,idx,itest,
3036	e t_seri, u_seri, v_seri,
3037	s zulow, zvlow, zustrdr, zvstrdr,
3038	s d_t_oro, d_u_oro, d_v_oro)
3039	ENDIF
3040	c
3041	c ajout des tendances
3042	!-----------------------------------------------------------------------------------------
3043	! ajout des tendances de la trainee de l'orographie
3044	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3045	!-----------------------------------------------------------------------------------------
3046	c
3047	ENDIF ! fin de test sur ok_orodr
3048	c
3049	if (mydebug) then
3050	call writefield_phy('u_seri',u_seri,llm)
3051	call writefield_phy('v_seri',v_seri,llm)
3052	call writefield_phy('t_seri',t_seri,llm)
3053	call writefield_phy('q_seri',q_seri,llm)
3054	endif
3055
3056	IF (ok_orolf) THEN
3057	c
3058	c selection des points pour lesquels le shema est actif:
3059	igwd=0
3060	DO i=1,klon
3061	itest(i)=0
3062	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3063	itest(i)=1
3064	igwd=igwd+1
3065	idx(igwd)=i
3066	ENDIF
3067	ENDDO
3068	c igwdim=MAX(1,igwd)
3069	c
3070	IF (ok_strato) THEN
3071
3072	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3073	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3074	e igwd,idx,itest,
3075	e t_seri, u_seri, v_seri,
3076	s zulow, zvlow, zustrli, zvstrli,
3077	s d_t_lif, d_u_lif, d_v_lif )
3078
3079	ELSE
3080	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3081	e rlat,zmea,zstd,zpic,
3082	e itest,
3083	e t_seri, u_seri, v_seri,
3084	s zulow, zvlow, zustrli, zvstrli,
3085	s d_t_lif, d_u_lif, d_v_lif)
3086	ENDIF
3087	c
3088	!-----------------------------------------------------------------------------------------
3089	! ajout des tendances de la portance de l'orographie
3090	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3091	!-----------------------------------------------------------------------------------------
3092	c
3093	ENDIF ! fin de test sur ok_orolf
3094	C HINES GWD PARAMETRIZATION
3095
3096	IF (ok_hines) then
3097
3098	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3099	i rlat,t_seri,u_seri,v_seri,
3100	o zustrhi,zvstrhi,
3101	o d_t_hin, d_u_hin, d_v_hin)
3102	c
3103	c ajout des tendances
3104	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3105
3106	ENDIF
3107	c
3108
3109	c
3110	cIM cf. FLott BEG
3111	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3112
3113	if (mydebug) then
3114	call writefield_phy('u_seri',u_seri,llm)
3115	call writefield_phy('v_seri',v_seri,llm)
3116	call writefield_phy('t_seri',t_seri,llm)
3117	call writefield_phy('q_seri',q_seri,llm)
3118	endif
3119
3120	DO i = 1, klon
3121	zustrph(i)=0.
3122	zvstrph(i)=0.
3123	ENDDO
3124	DO k = 1, klev
3125	DO i = 1, klon
3126	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3127	c (paprs(i,k)-paprs(i,k+1))/rg
3128	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3129	c (paprs(i,k)-paprs(i,k+1))/rg
3130	ENDDO
3131	ENDDO
3132	c
3133	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3134	c
3135	IF (is_sequential) THEN
3136
3137	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3138	C ra,rg,romega,
3139	C rlat,rlon,pphis,
3140	C zustrdr,zustrli,zustrph,
3141	C zvstrdr,zvstrli,zvstrph,
3142	C paprs,u,v,
3143	C aam, torsfc)
3144	ENDIF
3145	cIM cf. FLott END
3146	cIM
3147	IF (ip_ebil_phy.ge.2) THEN
3148	ztit='after orography'
3149	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3150	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3151	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3152	END IF
3153	c
3154	c
3155	cAA
3156	cAA Installation de l'interface online-offline pour traceurs
3157	cAA
3158	c====================================================================
3159	c Calcul des tendances traceurs
3160	c====================================================================
3161	C
3162
3163	call phytrac (
3164	I itap, days_elapsed+1, jH_cur, debut,
3165	I lafin, dtime, u, v, t,
3166	I paprs, pplay, pmfu, pmfd,
3167	I pen_u, pde_u, pen_d, pde_d,
3168	I cdragh, coefh, fm_therm, entr_therm,
3169	I u1, v1, ftsol, pctsrf,
3170	I rlat, frac_impa, frac_nucl,rlon,
3171	I presnivs, pphis, pphi, albsol1,
3172	I qx(:,:,ivap),rhcl, cldfra, rneb,
3173	I diafra, cldliq, itop_con, ibas_con,
3174	I pmflxr, pmflxs, prfl, psfl,
3175	I da, phi, mp, upwd,
3176	I dnwd, aerosol_couple, flxmass_w,
3177	I tau_aero, piz_aero, cg_aero, ccm,
3178	I rfname,
3179	O tr_seri)
3180
3181	IF (offline) THEN
3182
3183	print*,'Attention on met a 0 les thermiques pour phystoke'
3184	call phystokenc (
3185	I nlon,klev,pdtphys,rlon,rlat,
3186	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3187	I fm_therm,entr_therm,
3188	I cdragh,coefh,u1,v1,ftsol,pctsrf,
3189	I frac_impa, frac_nucl,
3190	I pphis,airephy,dtime,itap)
3191
3192
3193	ENDIF
3194
3195	c
3196	c Calculer le transport de l'eau et de l'energie (diagnostique)
3197	c
3198	CALL transp (paprs,zxtsol,
3199	e t_seri, q_seri, u_seri, v_seri, zphi,
3200	s ve, vq, ue, uq)
3201	c
3202	cIM global posePB BEG
3203	IF(1.EQ.0) THEN
3204	c
3205	CALL transp_lay (paprs,zxtsol,
3206	e t_seri, q_seri, u_seri, v_seri, zphi,
3207	s ve_lay, vq_lay, ue_lay, uq_lay)
3208	c
3209	ENDIF !(1.EQ.0) THEN
3210	cIM global posePB END
3211	c Accumuler les variables a stocker dans les fichiers histoire:
3212	c
3213	c+jld ec_conser
3214	DO k = 1, klev
3215	DO i = 1, klon
3216	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3217	d_t_ec(i,k)=0.5/ZRCPD
3218	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3219	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3220	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3221	END DO
3222	END DO
3223	c-jld ec_conser
3224	cIM
3225	IF (ip_ebil_phy.ge.1) THEN
3226	ztit='after physic'
3227	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3228	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3229	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3230	C Comme les tendances de la physique sont ajoute dans la dynamique,
3231	C on devrait avoir que la variation d'entalpie par la dynamique
3232	C est egale a la variation de la physique au pas de temps precedent.
3233	C Donc la somme de ces 2 variations devrait etre nulle.
3234	call diagphy(airephy,ztit,ip_ebil_phy
3235	e , topsw, toplw, solsw, sollw, sens
3236	e , evap, rain_fall, snow_fall, ztsol
3237	e , d_h_vcol, d_qt, d_ec
3238	s , fs_bound, fq_bound )
3239	C
3240	d_h_vcol_phy=d_h_vcol
3241	C
3242	END IF
3243	C
3244	c=======================================================================
3245	c SORTIES
3246	c=======================================================================
3247
3248	cIM Interpolation sur les niveaux de pression du NMC
3249	c -------------------------------------------------
3250	c
3251	#include "calcul_STDlev.h"
3252	twriteSTD(:,:,1)=tsumSTD(:,:,2)
3253	qwriteSTD(:,:,1)=qsumSTD(:,:,2)
3254	rhwriteSTD(:,:,1)=rhsumSTD(:,:,2)
3255	phiwriteSTD(:,:,1)=phisumSTD(:,:,2)
3256	uwriteSTD(:,:,1)=usumSTD(:,:,2)
3257	vwriteSTD(:,:,1)=vsumSTD(:,:,2)
3258	wwriteSTD(:,:,1)=wsumSTD(:,:,2)
3259
3260	twriteSTD(:,:,2)=tsumSTD(:,:,1)
3261	qwriteSTD(:,:,2)=qsumSTD(:,:,1)
3262	rhwriteSTD(:,:,2)=rhsumSTD(:,:,1)
3263	phiwriteSTD(:,:,2)=phisumSTD(:,:,1)
3264	uwriteSTD(:,:,2)=usumSTD(:,:,1)
3265	vwriteSTD(:,:,2)=vsumSTD(:,:,1)
3266	wwriteSTD(:,:,2)=wsumSTD(:,:,1)
3267
3268	twriteSTD(:,:,3)=tlevSTD(:,:)
3269	qwriteSTD(:,:,3)=qlevSTD(:,:)
3270	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3271	phiwriteSTD(:,:,3)=philevSTD(:,:)
3272	uwriteSTD(:,:,3)=ulevSTD(:,:)
3273	vwriteSTD(:,:,3)=vlevSTD(:,:)
3274	wwriteSTD(:,:,3)=wlevSTD(:,:)
3275
3276	twriteSTD(:,:,4)=tlevSTD(:,:)
3277	qwriteSTD(:,:,4)=qlevSTD(:,:)
3278	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3279	phiwriteSTD(:,:,4)=philevSTD(:,:)
3280	uwriteSTD(:,:,4)=ulevSTD(:,:)
3281	vwriteSTD(:,:,4)=vlevSTD(:,:)
3282	wwriteSTD(:,:,4)=wlevSTD(:,:)
3283	c
3284	c slp sea level pressure
3285	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3286	c
3287	ccc prw = eau precipitable
3288	DO i = 1, klon
3289	prw(i) = 0.
3290	DO k = 1, klev
3291	prw(i) = prw(i) +
3292	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3293	ENDDO
3294	ENDDO
3295	c
3296	cIM initialisation + calculs divers diag AMIP2
3297	c
3298	#include "calcul_divers.h"
3299	c
3300	IF (config_inca /= 'none') THEN
3301	#ifdef INCA
3302	CALL VTe(VTphysiq)
3303	CALL VTb(VTinca)
3304
3305	CALL chemhook_end (calday,
3306	$ dtime,
3307	$ pplay,
3308	$ t_seri,
3309	$ tr_seri,
3310	$ nbtr,
3311	$ paprs,
3312	$ q_seri,
3313	$ annee_ref,
3314	$ day_ini,
3315	$ airephy,
3316	$ xjour,
3317	$ pphi,
3318	$ pphis,
3319	$ zx_rh)
3320
3321	CALL VTe(VTinca)
3322	CALL VTb(VTphysiq)
3323	#endif
3324	END IF
3325
3326	c=============================================================
3327	c
3328	c Convertir les incrementations en tendances
3329	c
3330	if (mydebug) then
3331	call writefield_phy('u_seri',u_seri,llm)
3332	call writefield_phy('v_seri',v_seri,llm)
3333	call writefield_phy('t_seri',t_seri,llm)
3334	call writefield_phy('q_seri',q_seri,llm)
3335	endif
3336
3337	DO k = 1, klev
3338	DO i = 1, klon
3339	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3340	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3341	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3342	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3343	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3344	ENDDO
3345	ENDDO
3346	c
3347	IF (nqtot.GE.3) THEN
3348	DO iq = 3, nqtot
3349	DO k = 1, klev
3350	DO i = 1, klon
3351	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3352	ENDDO
3353	ENDDO
3354	ENDDO
3355	ENDIF
3356	c
3357	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3358	cIM global posePB#include "write_bilKP_ins.h"
3359	cIM global posePB#include "write_bilKP_ave.h"
3360	c
3361	c Sauvegarder les valeurs de t et q a la fin de la physique:
3362	c
3363	DO k = 1, klev
3364	DO i = 1, klon
3365	u_ancien(i,k) = u_seri(i,k)
3366	v_ancien(i,k) = v_seri(i,k)
3367	t_ancien(i,k) = t_seri(i,k)
3368	q_ancien(i,k) = q_seri(i,k)
3369	ENDDO
3370	ENDDO
3371	c
3372	!==========================================================================
3373	! Sorties des tendances pour un point particulier
3374	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3375	! pour le debug
3376	! La valeur de igout est attribuee plus haut dans le programme
3377	!==========================================================================
3378
3379	if (prt_level.ge.1) then
3380	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3381	write(lunout,*)
3382	s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3383	write(lunout,*)
3384	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
3385	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3386	s pctsrf(igout,is_sic)
3387	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3388	do k=1,klev
3389	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3390	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3391	s d_t_eva(igout,k)
3392	enddo
3393	write(lunout,*) 'cool,heat'
3394	do k=1,klev
3395	write(lunout,*) cool(igout,k),heat(igout,k)
3396	enddo
3397
3398	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3399	do k=1,klev
3400	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3401	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3402	enddo
3403
3404	write(lunout,*) 'd_ps ',d_ps(igout)
3405	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3406	do k=1,klev
3407	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3408	s d_qx(igout,k,1),d_qx(igout,k,2)
3409	enddo
3410	endif
3411
3412	!==========================================================================
3413
3414	c============================================================
3415	c Calcul de la temperature potentielle
3416	c============================================================
3417	DO k = 1, klev
3418	DO i = 1, klon
3419	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3420	ENDDO
3421	ENDDO
3422	c
3423
3424	c 22.03.04 BEG
3425	c=============================================================
3426	c Ecriture des sorties
3427	c=============================================================
3428	#ifdef CPP_IOIPSL
3429
3430	c Recupere des varibles calcule dans differents modules
3431	c pour ecriture dans histxxx.nc
3432
3433	! Get some variables from module fonte_neige_mod
3434	CALL fonte_neige_get_vars(pctsrf,
3435	. zxfqcalving, zxfqfonte, zxffonte)
3436
3437
3438	#include "phys_output_write.h"
3439
3440	#ifdef histISCCP
3441	#include "write_histISCCP.h"
3442	#endif
3443
3444	#ifdef histmthNMC
3445	#include "write_histmthNMC.h"
3446	#endif
3447
3448	#include "write_histday_seri.h"
3449
3450	#include "write_paramLMDZ_phy.h"
3451
3452	#endif
3453
3454	c 22.03.04 END
3455	c
3456	c====================================================================
3457	c Si c'est la fin, il faut conserver l'etat de redemarrage
3458	c====================================================================
3459	c
3460
3461
3462	IF (lafin) THEN
3463	itau_phy = itau_phy + itap
3464	CALL phyredem ("restartphy.nc")
3465	! open(97,form="unformatted",file="finbin")
3466	! write(97) u_seri,v_seri,t_seri,q_seri
3467	! close(97)
3468	C$OMP single
3469	if (read_climoz) then
3470	if (is_mpi_root) then
3471	call nf95_close(ncid_climoz)
3472	end if
3473	deallocate(press_climoz) ! pointer
3474	end if
3475	C$OMP end single nowait
3476	ENDIF
3477
3478	! first=.false.
3479
3480	RETURN
3481	END
3482	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3483	IMPLICIT none
3484	c
3485	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3486	c la conservation de l'eau
3487	c
3488	#include "YOMCST.h"
3489	INTEGER klon,klev
3490	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3491	REAL aire(klon)
3492	REAL qtotal, zx, qcheck
3493	INTEGER i, k
3494	c
3495	zx = 0.0
3496	DO i = 1, klon
3497	zx = zx + aire(i)
3498	ENDDO
3499	qtotal = 0.0
3500	DO k = 1, klev
3501	DO i = 1, klon
3502	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3503	. *(paprs(i,k)-paprs(i,k+1))/RG
3504	ENDDO
3505	ENDDO
3506	c
3507	qcheck = qtotal/zx
3508	c
3509	RETURN
3510	END
3511	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3512	IMPLICIT none
3513	c
3514	c Tranformer une variable de la grille physique a
3515	c la grille d'ecriture
3516	c
3517	INTEGER nfield,nlon,iim,jjmp1, jjm
3518	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3519	c
3520	INTEGER i, n, ig
3521	c
3522	jjm = jjmp1 - 1
3523	DO n = 1, nfield
3524	DO i=1,iim
3525	ecrit(i,n) = fi(1,n)
3526	ecrit(i+jjm*iim,n) = fi(nlon,n)
3527	ENDDO
3528	DO ig = 1, nlon - 2
3529	ecrit(iim+ig,n) = fi(1+ig,n)
3530	ENDDO
3531	ENDDO
3532	RETURN
3533	END

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

Download in other formats:

Original Format