Context Navigation

physiq.F @ 1391

Last change on this file since 1391 was 1391, checked in by musat, 14 years ago
Variables (read_climoz, ncid_climoz, press_climoz) are "OpenMP shared" and do not need $OMP THREADPRIVATE. LG/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 120.2 KB

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

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

Download in other formats:

Original Format