Context Navigation

physiq.F @ 1244

Last change on this file since 1244 was 1244, checked in by Laurent Fairhead, 15 years ago
Le calcul de ratqs faisait osciller le modele et donnait une solution asymptotique bidon et dépendant fortement du pas de temps.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 115.1 KB

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

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

Download in other formats:

Original Format