Context Navigation

physiq.F @ 1322

Last change on this file since 1322 was 1322, checked in by Laurent Fairhead, 14 years ago

Improvements concerning wake parametrisation (from JYG, NR, IT, with more to come).
Alp_offset is read in form physiq.def file

Améliorations à la paramétrisation des poches froides (de JYG, NR, IT, d'autres
sont à venir)
Alp_offset est rajouté à la liste des paramètres lus dans physiq.def

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

File size: 119.8 KB

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

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

Download in other formats:

Original Format