Context Navigation

physiq.F @ 1215

Last change on this file since 1215 was 1215, checked in by lguez, 15 years ago

-- Made "ozonecm" a function instead of a subroutine. Used assumed shape
arguments in "ozonecm".

-- Corrected long name and computation of NetCDF variable "ozone" in
the files "hist*".

-- Corrected comments for ozone variables.

-- In the case "read_climoz", used variables "rmd" and "rmo3" from
"YOMCST.h" instead of writing approximate values.

-- Replaced "real*..." declarations (not conforming to Fortran standard)
by "real(kind=...)" declarations.

-- Replaced value "1./46.6968" in ozone computations by the equivalent
(but clearer) "dobson_u * 1e3" (relative difference ~ 1e-5).

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

File size: 114.1 KB

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

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

Download in other formats:

Original Format