Context Navigation

physiq.F @ 1191

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

Reecriture de phytrac et les routines concernes (Anthony Jamelot)

les suffix change de F -> F90 (nflxtr.F90,cltracrn.F90,initrrnpb.F90,cvltr.F90,minmaxqfi.F90,cltrac.F90,phytrac.F90)

Traitement d'un nouveau traceur berelium (optionel, toujours pour des
tests)(Anthony Jamelot)

radiornpb.F change du nom pour radio_decay.F90 car il traite maintenant tout les traceurs radioactives
ajoute init_be.F90

Nouveau interface dans phytrac pour serparer les calculs et appels
specifique a INCA avec les traitements des traceurs specifiques au LMDZ
(JG)

ajoute tracinca_mod.F90 pour les appeles a INCA
ajoute traclmdz_mod.F90 pour les calculs des traceurs specifiques a LMDZ
enleve fichier restartrac et ajoute la variable trs dans restartphy.nc

La convergence numerique a etait rompue uniquement pour les traceurs
LMDZ RN et PB.

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

File size: 112.9 KB

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

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

Download in other formats:

Original Format