Context Navigation

physiq.F @ 1179

Last change on this file since 1179 was 1179, checked in by jghattas, 15 years ago

Ajout de l'interpolation vertical pour les nouveaux fichiers de forcage des aerosols. Utilisant les anciennes fichiers de SO4 pas d'interpolation possible. Convergence numerique avec la version precedente en utilisant les anciens fichiers des SO4. aerosol_optic.F90 change du nom pour readaerosol_optic.F90 (lecture d'aerosol + optic) Les fichiers de forcage aerosol doit maintenant avoir le suffix .nc.

Correction des bugs pour inca et certain diagnostiques optionelles de radlwsw.
Ajout de test pour le choix advection schema.

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

File size: 114.2 KB

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

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

Download in other formats:

Original Format