Context Navigation

physiq.F @ 1355

Last change on this file since 1355 was 1355, checked in by musat, 14 years ago
use of phys_cal_mod in conf_phys (and change call order in physiq.F) to automatically calculate the output frequency of standard pressure monthly file (histmthNMC.nc) 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 1355 2010-04-14 13:27:18Z 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	cIM for NMC files en parallele
1130	c$OMP THREADPRIVATE(read_climoz, ncid_climoz, press_climoz)
1131
1132	integer, save:: co3i = 0
1133	! time index in NetCDF file of current ozone fields
1134	c$OMP THREADPRIVATE(co3i)
1135
1136	integer ro3i
1137	! required time index in NetCDF file for the ozone fields, between 1
1138	! and 360
1139
1140	#include "YOMCST.h"
1141	#include "YOETHF.h"
1142	#include "FCTTRE.h"
1143	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1144	#include "conema3.h"
1145	#include "fisrtilp.h"
1146	#include "nuage.h"
1147	#include "compbl.h"
1148	cIM 100106 END : pouvoir sortir les ctes de la physique
1149	c
1150	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1151	c Declarations pour Simulateur COSP
1152	c============================================================
1153	real :: mr_ozone(klon,klev)
1154	cIM for NMC files
1155	missing_val=nf90_fill_real
1156	c======================================================================
1157	! Gestion calendrier : mise a jour du module phys_cal_mod
1158	!
1159	CALL phys_cal_update(jD_cur,jH_cur)
1160
1161	c======================================================================
1162	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1163	! Utilise notamment en 1D mais peut etre active egalement en 3D
1164	! en imposant la valeur de igout.
1165	c======================================================================
1166
1167	if (prt_level.ge.1) then
1168	igout=klon/2+1/klon
1169	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1170	write(lunout,*)
1171	s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1172	write(lunout,*)
1173	s nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1174
1175	write(lunout,*) 'paprs, play, phi, u, v, t'
1176	do k=1,klev
1177	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1178	s u(igout,k),v(igout,k),t(igout,k)
1179	enddo
1180	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1181	do k=1,klev
1182	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1183	enddo
1184	endif
1185
1186	c======================================================================
1187
1188	cym => necessaire pour iflag_con != 2
1189	pmfd(:,:) = 0.
1190	pen_u(:,:) = 0.
1191	pen_d(:,:) = 0.
1192	pde_d(:,:) = 0.
1193	pde_u(:,:) = 0.
1194	aam=0.
1195
1196	torsfc=0.
1197	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1198
1199	if (first) then
1200
1201	cCR:nvelles variables convection/poches froides
1202
1203	print*, '================================================='
1204	print*, 'Allocation des variables locales et sauvegardees'
1205	call phys_local_var_init
1206	c
1207	pasphys=pdtphys
1208	c appel a la lecture du run.def physique
1209	call conf_phys(ok_journe, ok_mensuel,
1210	. ok_instan, ok_hf,
1211	. ok_LES,
1212	. solarlong0,seuil_inversion,
1213	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1214	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
1215	. ok_ade, ok_aie, aerosol_couple,
1216	. flag_aerosol, new_aod,
1217	. bl95_b0, bl95_b1,
1218	. iflag_thermals,nsplit_thermals,tau_thermals,
1219	. iflag_thermals_ed,iflag_thermals_optflux,
1220	c nv flags pour la convection et les poches froides
1221	. iflag_coupl,iflag_clos,iflag_wake, read_climoz)
1222	call phys_state_var_init(read_climoz)
1223	call phys_output_var_init
1224	print*, '================================================='
1225	cIM for NMC files
1226	cIM freq_moyNMC = frequences auxquelles on moyenne les champs accumules
1227	cIM sur les niveaux de pression standard du NMC
1228	DO n=1, nout
1229	freq_moyNMC(n)=freq_outNMC(n)/freq_calNMC(n)
1230	ENDDO
1231	c
1232	cIM beg
1233	dnwd0=0.0
1234	ftd=0.0
1235	fqd=0.0
1236	cin=0.
1237	cym Attention pbase pas initialise dans concvl !!!!
1238	pbase=0
1239	paire_ter(:)=0.
1240	cIM 180608
1241	c pmflxr=0.
1242	c pmflxs=0.
1243	itau_con=0
1244	first=.false.
1245
1246	endif ! first
1247
1248	modname = 'physiq'
1249	cIM
1250	IF (ip_ebil_phy.ge.1) THEN
1251	DO i=1,klon
1252	zero_v(i)=0.
1253	END DO
1254	END IF
1255	ok_sync=.TRUE.
1256
1257	IF (debut) THEN
1258	CALL suphel ! initialiser constantes et parametres phys.
1259	ENDIF
1260
1261	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1262
1263
1264	c======================================================================
1265	! Gestion calendrier : mise a jour du module phys_cal_mod
1266	!
1267	cIM CALL phys_cal_update(jD_cur,jH_cur)
1268
1269	c
1270	c Si c'est le debut, il faut initialiser plusieurs choses
1271	c ********
1272	c
1273	IF (debut) THEN
1274	!rv
1275	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1276	cde la convection a partir des caracteristiques du thermique
1277	wght_th(:,:)=1.
1278	lalim_conv(:)=1
1279	cRC
1280	u10m(:,:)=0.
1281	v10m(:,:)=0.
1282	rain_con(:)=0.
1283	snow_con(:)=0.
1284	topswai(:)=0.
1285	topswad(:)=0.
1286	solswai(:)=0.
1287	solswad(:)=0.
1288
1289	lambda_th(:,:)=0.
1290	wmax_th(:)=0.
1291	tau_overturning_th(:)=0.
1292
1293	IF (config_inca /= 'none') THEN
1294	! jg : initialisation jusqu'au ces variables sont dans restart
1295	ccm(:,:,:) = 0.
1296	tau_aero(:,:,:,:) = 0.
1297	piz_aero(:,:,:,:) = 0.
1298	cg_aero(:,:,:,:) = 0.
1299	END IF
1300
1301	rnebcon0(:,:) = 0.0
1302	clwcon0(:,:) = 0.0
1303	rnebcon(:,:) = 0.0
1304	clwcon(:,:) = 0.0
1305
1306	cIM
1307	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1308	c
1309	print*,'iflag_coupl,iflag_clos,iflag_wake',
1310	. iflag_coupl,iflag_clos,iflag_wake
1311	print*,'CYCLE_DIURNE', cycle_diurne
1312	c
1313	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1314	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1315	CALL abort_gcm (modname,abort_message,1)
1316	ENDIF
1317	c
1318	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1319	abort_message = 'ISCCP-like outputs may be available for KE
1320	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1321	CALL abort_gcm (modname,abort_message,1)
1322	ENDIF
1323	c
1324	c Initialiser les compteurs:
1325	c
1326	itap = 0
1327	itaprad = 0
1328
1329	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1330	!! Un petit travail à faire ici.
1331	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1332
1333	if (iflag_pbl>1) then
1334	PRINT*, "Using method MELLOR&YAMADA"
1335	endif
1336
1337	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1338	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1339	! dyn3d
1340	! Attention : la version precedente n'etait pas tres propre.
1341	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1342	! pour obtenir le meme resultat.
1343	dtime=pdtphys
1344	radpas = NINT( 86400./dtime/nbapp_rad)
1345	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1346
1347	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1348	cIM begin
1349	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1350	$,ratqs(1,1)
1351	cIM end
1352
1353
1354
1355	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1356	c
1357	C on remet le calendrier a zero
1358	c
1359	IF (raz_date .eq. 1) THEN
1360	itau_phy = 0
1361	ENDIF
1362
1363	cIM cf. AM 081204 BEG
1364	PRINT*,'cycle_diurne3 =',cycle_diurne
1365	cIM cf. AM 081204 END
1366	c
1367	CALL printflag( tabcntr0,radpas,ok_journe,
1368	, ok_instan, ok_region )
1369	c
1370	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1371	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1372	. pdtphys
1373	abort_message='Pas physique n est pas correct '
1374	! call abort_gcm(modname,abort_message,1)
1375	dtime=pdtphys
1376	ENDIF
1377	IF (nlon .NE. klon) THEN
1378	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1379	. klon
1380	abort_message='nlon et klon ne sont pas coherents'
1381	call abort_gcm(modname,abort_message,1)
1382	ENDIF
1383	IF (nlev .NE. klev) THEN
1384	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1385	. klev
1386	abort_message='nlev et klev ne sont pas coherents'
1387	call abort_gcm(modname,abort_message,1)
1388	ENDIF
1389	c
1390	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1391	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1392	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1393	abort_message='Nbre d appels au rayonnement insuffisant'
1394	call abort_gcm(modname,abort_message,1)
1395	ENDIF
1396	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1397	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1398	. ok_cvl
1399	c
1400	cKE43
1401	c Initialisation pour la convection de K.E. (sb):
1402	IF (iflag_con.GE.3) THEN
1403
1404	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1405	WRITE(lunout,*)
1406	. "On va utiliser le melange convectif des traceurs qui"
1407	WRITE(lunout,*)"est calcule dans convect4.3"
1408	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1409
1410	DO i = 1, klon
1411	ema_cbmf(i) = 0.
1412	ema_pcb(i) = 0.
1413	ema_pct(i) = 0.
1414	ema_workcbmf(i) = 0.
1415	ENDDO
1416	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1417	DO i = 1, klon
1418	ibas_con(i) = 1
1419	itop_con(i) = 1
1420	ENDDO
1421	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1422	c===============================================================================
1423	cCR:04.12.07: initialisations poches froides
1424	c Controle de ALE et ALP pour la fermeture convective (jyg)
1425	if (iflag_wake.eq.1) then
1426	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1427	s ,alp_bl_prescr, ale_bl_prescr)
1428	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1429	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1430	endif
1431
1432	do i = 1,klon
1433	Ale_bl(i)=0.
1434	Alp_bl(i)=0.
1435	enddo
1436
1437	c================================================================================
1438
1439	ENDIF
1440
1441	DO i=1,klon
1442	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1443	ENDDO
1444
1445	c34EK
1446	IF (ok_orodr) THEN
1447
1448	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1449	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1450	! justement quand ok_orodr = false.
1451	! ce rugoro est utilise par la couche limite et fait double emploi
1452	! avec les paramétrisations spécifiques de Francois Lott.
1453	! DO i=1,klon
1454	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1455	! ENDDO
1456	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1457	IF (ok_strato) THEN
1458	CALL SUGWD_strato(klon,klev,paprs,pplay)
1459	ELSE
1460	CALL SUGWD(klon,klev,paprs,pplay)
1461	ENDIF
1462
1463	DO i=1,klon
1464	zuthe(i)=0.
1465	zvthe(i)=0.
1466	if(zstd(i).gt.10.)then
1467	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1468	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1469	endif
1470	ENDDO
1471	ENDIF
1472	c
1473	c
1474	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1475	WRITE(lunout,*)'La frequence de lecture surface est de ',
1476	. lmt_pas
1477	c
1478	cIM 030306 END
1479
1480	capemaxcels = 't_max(X)'
1481	t2mincels = 't_min(X)'
1482	t2maxcels = 't_max(X)'
1483	tinst = 'inst(X)'
1484	tave = 'ave(X)'
1485	cIM cf. AM 081204 BEG
1486	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1487	cIM cf. AM 081204 END
1488	c
1489	c=============================================================
1490	c Initialisation des sorties
1491	c=============================================================
1492
1493	#ifdef CPP_IOIPSL
1494
1495	c$OMP MASTER
1496	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1497	& ctetaSTD,dtime,ok_veget,
1498	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1499	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie,
1500	& read_climoz, new_aod, aerosol_couple)
1501	c$OMP END MASTER
1502	c$OMP BARRIER
1503
1504	#ifdef histISCCP
1505	#include "ini_histISCCP.h"
1506	#endif
1507
1508	#ifdef histNMC
1509	#include "ini_histhfNMC.h"
1510	#include "ini_histdayNMC.h"
1511	#include "ini_histmthNMC.h"
1512	#endif
1513
1514	#include "ini_histday_seri.h"
1515
1516	#include "ini_paramLMDZ_phy.h"
1517
1518	#endif
1519
1520	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1521	ecrit_hf2mth = ecrit_mth/ecrit_hf
1522
1523	ecrit_hf = ecrit_hf * un_jour
1524	cIM
1525	IF(ecrit_day.LE.1.) THEN
1526	ecrit_day = ecrit_day * un_jour !en secondes
1527	ENDIF
1528	cIM
1529	ecrit_mth = ecrit_mth * un_jour
1530	ecrit_ins = ecrit_ins * un_jour
1531	ecrit_reg = ecrit_reg * un_jour
1532	ecrit_tra = ecrit_tra * un_jour
1533	ecrit_ISCCP = ecrit_ISCCP * un_jour
1534	ecrit_LES = ecrit_LES * un_jour
1535	c
1536	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1537	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1538	. ecrit_hf2mth
1539	cIM 030306 END
1540
1541
1542	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1543	date0 = jD_ref
1544	WRITE(,) 'physiq date0 : ',date0
1545	c
1546	c
1547	c
1548	c Prescrire l'ozone dans l'atmosphere
1549	c
1550	c
1551	cc DO i = 1, klon
1552	cc DO k = 1, klev
1553	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1554	cc ENDDO
1555	cc ENDDO
1556	c
1557	IF (config_inca /= 'none') THEN
1558	#ifdef INCA
1559	CALL VTe(VTphysiq)
1560	CALL VTb(VTinca)
1561	! iii = MOD(NINT(xjour),360)
1562	! calday = FLOAT(iii) + jH_cur
1563	calday = FLOAT(days_elapsed) + jH_cur
1564	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1565
1566	CALL chemini(
1567	$ rg,
1568	$ ra,
1569	$ airephy,
1570	$ rlat,
1571	$ rlon,
1572	$ presnivs,
1573	$ calday,
1574	$ klon,
1575	$ nqtot,
1576	$ pdtphys,
1577	$ annee_ref,
1578	$ day_ref,
1579	$ itau_phy)
1580
1581	CALL VTe(VTinca)
1582	CALL VTb(VTphysiq)
1583	#endif
1584	END IF
1585	c
1586	c
1587	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1588	! Nouvelle initialisation pour le rayonnement RRTM
1589	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1590
1591	call iniradia(klon,klev,paprs(1,1:klev+1))
1592
1593	C$omp single
1594	if (read_climoz >= 1) then
1595	call open_climoz(ncid_climoz, press_climoz)
1596	END IF
1597	C$omp end single
1598	ENDIF
1599	!
1600	! ************** Fin de IF ( debut ) *************
1601	!
1602	!
1603	! Incrementer le compteur de la physique
1604	!
1605	itap = itap + 1
1606	!
1607	! Update fraction of the sub-surfaces (pctsrf) and
1608	! initialize, where a new fraction has appeared, all variables depending
1609	! on the surface fraction.
1610	!
1611	CALL change_srf_frac(itap, dtime, days_elapsed+1,
1612	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1613
1614	! Tendances bidons pour les processus qui n'affectent pas certaines
1615	! variables.
1616	du0(:,:)=0.
1617	dv0(:,:)=0.
1618	dq0(:,:)=0.
1619	dql0(:,:)=0.
1620	c
1621	c Mettre a zero des variables de sortie (pour securite)
1622	c
1623	DO i = 1, klon
1624	d_ps(i) = 0.0
1625	ENDDO
1626	DO k = 1, klev
1627	DO i = 1, klon
1628	d_t(i,k) = 0.0
1629	d_u(i,k) = 0.0
1630	d_v(i,k) = 0.0
1631	ENDDO
1632	ENDDO
1633	DO iq = 1, nqtot
1634	DO k = 1, klev
1635	DO i = 1, klon
1636	d_qx(i,k,iq) = 0.0
1637	ENDDO
1638	ENDDO
1639	ENDDO
1640	da(:,:)=0.
1641	mp(:,:)=0.
1642	phi(:,:,:)=0.
1643	c
1644	c Ne pas affecter les valeurs entrees de u, v, h, et q
1645	c
1646	DO k = 1, klev
1647	DO i = 1, klon
1648	t_seri(i,k) = t(i,k)
1649	u_seri(i,k) = u(i,k)
1650	v_seri(i,k) = v(i,k)
1651	q_seri(i,k) = qx(i,k,ivap)
1652	ql_seri(i,k) = qx(i,k,iliq)
1653	qs_seri(i,k) = 0.
1654	ENDDO
1655	ENDDO
1656	IF (nqtot.GE.3) THEN
1657	DO iq = 3, nqtot
1658	DO k = 1, klev
1659	DO i = 1, klon
1660	tr_seri(i,k,iq-2) = qx(i,k,iq)
1661	ENDDO
1662	ENDDO
1663	ENDDO
1664	ELSE
1665	DO k = 1, klev
1666	DO i = 1, klon
1667	tr_seri(i,k,1) = 0.0
1668	ENDDO
1669	ENDDO
1670	ENDIF
1671	C
1672	DO i = 1, klon
1673	ztsol(i) = 0.
1674	ENDDO
1675	DO nsrf = 1, nbsrf
1676	DO i = 1, klon
1677	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1678	ENDDO
1679	ENDDO
1680	cIM
1681	IF (ip_ebil_phy.ge.1) THEN
1682	ztit='after dynamic'
1683	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1684	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1685	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1686	C Comme les tendances de la physique sont ajoute dans la dynamique,
1687	C on devrait avoir que la variation d'entalpie par la dynamique
1688	C est egale a la variation de la physique au pas de temps precedent.
1689	C Donc la somme de ces 2 variations devrait etre nulle.
1690	call diagphy(airephy,ztit,ip_ebil_phy
1691	e , zero_v, zero_v, zero_v, zero_v, zero_v
1692	e , zero_v, zero_v, zero_v, ztsol
1693	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1694	s , fs_bound, fq_bound )
1695	END IF
1696
1697	c Diagnostiquer la tendance dynamique
1698	c
1699	IF (ancien_ok) THEN
1700	DO k = 1, klev
1701	DO i = 1, klon
1702	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1703	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1704	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1705	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1706	ENDDO
1707	ENDDO
1708	ELSE
1709	DO k = 1, klev
1710	DO i = 1, klon
1711	d_u_dyn(i,k) = 0.0
1712	d_v_dyn(i,k) = 0.0
1713	d_t_dyn(i,k) = 0.0
1714	d_q_dyn(i,k) = 0.0
1715	ENDDO
1716	ENDDO
1717	ancien_ok = .TRUE.
1718	ENDIF
1719	c
1720	c Ajouter le geopotentiel du sol:
1721	c
1722	DO k = 1, klev
1723	DO i = 1, klon
1724	zphi(i,k) = pphi(i,k) + pphis(i)
1725	ENDDO
1726	ENDDO
1727	c
1728	c Verifier les temperatures
1729	c
1730	cIM BEG
1731	IF (check) THEN
1732	amn=MIN(ftsol(1,is_ter),1000.)
1733	amx=MAX(ftsol(1,is_ter),-1000.)
1734	DO i=2, klon
1735	amn=MIN(ftsol(i,is_ter),amn)
1736	amx=MAX(ftsol(i,is_ter),amx)
1737	ENDDO
1738	c
1739	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1740	ENDIF !(check) THEN
1741	cIM END
1742	c
1743	CALL hgardfou(t_seri,ftsol,'debutphy')
1744	c
1745	cIM BEG
1746	IF (check) THEN
1747	amn=MIN(ftsol(1,is_ter),1000.)
1748	amx=MAX(ftsol(1,is_ter),-1000.)
1749	DO i=2, klon
1750	amn=MIN(ftsol(i,is_ter),amn)
1751	amx=MAX(ftsol(i,is_ter),amx)
1752	ENDDO
1753	c
1754	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1755	ENDIF !(check) THEN
1756	cIM END
1757	c
1758	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1759	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1760	c
1761	if (read_climoz >= 1) then
1762	C Ozone from a file
1763	! Update required ozone index:
1764	ro3i = int((days_elapsed + jh_cur - jh_1jan)
1765	$ / ioget_year_len(year_cur) * 360.) + 1
1766	if (ro3i == 361) ro3i = 360
1767	C (This should never occur, except perhaps because of roundup
1768	C error. See documentation.)
1769	if (ro3i /= co3i) then
1770	C Update ozone field:
1771	if (read_climoz == 1) then
1772	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
1773	$ press_in_edg=press_climoz, paprs=paprs, v3=wo)
1774	else
1775	C read_climoz == 2
1776	call regr_pr_av(ncid_climoz,
1777	$ (/"tro3 ", "tro3_daylight"/),
1778	$ julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
1779	$ v3=wo)
1780	end if
1781	! Convert from mole fraction of ozone to column density of ozone in a
1782	! cell, in kDU:
1783	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
1784	$ * rmo3 / rmd * zmasse / dobson_u / 1e3
1785	C (By regridding ozone values for LMDZ only once every 360th of
1786	C year, we have already neglected the variation of pressure in one
1787	C 360th of year. So do not recompute "wo" at each time step even if
1788	C "zmasse" changes a little.)
1789	co3i = ro3i
1790	end if
1791	elseif (MOD(itap-1,lmt_pas) == 0) THEN
1792	C Once per day, update ozone from Royer:
1793	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
1794	ENDIF
1795	c
1796	c Re-evaporer l'eau liquide nuageuse
1797	c
1798	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1799	DO i = 1, klon
1800	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1801	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1802	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1803	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1804	zb = MAX(0.0,ql_seri(i,k))
1805	za = - MAX(0.0,ql_seri(i,k))
1806	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1807	t_seri(i,k) = t_seri(i,k) + za
1808	q_seri(i,k) = q_seri(i,k) + zb
1809	ql_seri(i,k) = 0.0
1810	d_t_eva(i,k) = za
1811	d_q_eva(i,k) = zb
1812	ENDDO
1813	ENDDO
1814	cIM
1815	IF (ip_ebil_phy.ge.2) THEN
1816	ztit='after reevap'
1817	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1818	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1819	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1820	call diagphy(airephy,ztit,ip_ebil_phy
1821	e , zero_v, zero_v, zero_v, zero_v, zero_v
1822	e , zero_v, zero_v, zero_v, ztsol
1823	e , d_h_vcol, d_qt, d_ec
1824	s , fs_bound, fq_bound )
1825	C
1826	END IF
1827
1828	c
1829	c=========================================================================
1830	! Calculs de l'orbite.
1831	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1832	! doit donc etre placé avant radlwsw et pbl_surface
1833
1834	! calcul selon la routine utilisee pour les planetes
1835	if (new_orbit) then
1836	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1837	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1838	! day_since_equinox = (jD_cur) - jD_eq
1839	call solarlong(day_since_equinox, zlongi, dist)
1840	else
1841	! calcul selon la routine utilisee pour l'AR4
1842	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1843	! solarlong0
1844	if (solarlong0<-999.) then
1845	CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
1846	else
1847	zlongi=solarlong0 ! longitude solaire vraie
1848	dist=1. ! distance au soleil / moyenne
1849	endif
1850	endif
1851	if(prt_level.ge.1) &
1852	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1853
1854	! Avec ou sans cycle diurne
1855	IF (cycle_diurne) THEN
1856	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1857	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1858	ELSE
1859	CALL angle(zlongi, rlat, fract, rmu0)
1860	ENDIF
1861
1862	if (mydebug) then
1863	call writefield_phy('u_seri',u_seri,llm)
1864	call writefield_phy('v_seri',v_seri,llm)
1865	call writefield_phy('t_seri',t_seri,llm)
1866	call writefield_phy('q_seri',q_seri,llm)
1867	endif
1868
1869	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1870	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1871	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1872	c turbulent du couche limit.
1873	c
1874	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1875	c ecriture des fihiers hist_XXXX.nc, ces sont :
1876	c qsol, zq2m, s_pblh, s_lcl,
1877	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1878	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1879	c zxrugs, zu10m, zv10m, fder,
1880	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1881	c frugs, agesno, fsollw, fsolsw,
1882	c d_ts, fevap, fluxlat, t2m,
1883	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1884	c
1885	c Certains ne sont pas utiliser du tout :
1886	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1887	c
1888
1889	CALL pbl_surface(
1890	e dtime, date0, itap, days_elapsed+1,
1891	e debut, lafin,
1892	e rlon, rlat, rugoro, rmu0,
1893	e rain_fall, snow_fall, solsw, sollw,
1894	e t_seri, q_seri, u_seri, v_seri,
1895	e pplay, paprs, pctsrf,
1896	+ ftsol, falb1, falb2, u10m, v10m,
1897	s sollwdown, cdragh, cdragm, u1, v1,
1898	s albsol1, albsol2, sens, evap,
1899	s zxtsol, zxfluxlat, zt2m, qsat2m,
1900	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1901	s coefh, slab_wfbils,
1902	d qsol, zq2m, s_pblh, s_lcl,
1903	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1904	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1905	d zxrugs, zu10m, zv10m, fder,
1906	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1907	d frugs, agesno, fsollw, fsolsw,
1908	d d_ts, fevap, fluxlat, t2m,
1909	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1910	- dsens, devap, zxsnow,
1911	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1912
1913
1914	!-----------------------------------------------------------------------------------------
1915	! ajout des tendances de la diffusion turbulente
1916	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1917	!-----------------------------------------------------------------------------------------
1918
1919	if (mydebug) then
1920	call writefield_phy('u_seri',u_seri,llm)
1921	call writefield_phy('v_seri',v_seri,llm)
1922	call writefield_phy('t_seri',t_seri,llm)
1923	call writefield_phy('q_seri',q_seri,llm)
1924	endif
1925
1926
1927	IF (ip_ebil_phy.ge.2) THEN
1928	ztit='after surface_main'
1929	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1930	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1931	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1932	call diagphy(airephy,ztit,ip_ebil_phy
1933	e , zero_v, zero_v, zero_v, zero_v, sens
1934	e , evap , zero_v, zero_v, ztsol
1935	e , d_h_vcol, d_qt, d_ec
1936	s , fs_bound, fq_bound )
1937	END IF
1938
1939	c =================================================================== c
1940	c Calcul de Qsat
1941
1942	DO k = 1, klev
1943	DO i = 1, klon
1944	zx_t = t_seri(i,k)
1945	IF (thermcep) THEN
1946	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1947	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1948	zx_qs = MIN(0.5,zx_qs)
1949	zcor = 1./(1.-retv*zx_qs)
1950	zx_qs = zx_qs*zcor
1951	ELSE
1952	IF (zx_t.LT.t_coup) THEN
1953	zx_qs = qsats(zx_t)/pplay(i,k)
1954	ELSE
1955	zx_qs = qsatl(zx_t)/pplay(i,k)
1956	ENDIF
1957	ENDIF
1958	zqsat(i,k)=zx_qs
1959	ENDDO
1960	ENDDO
1961
1962	if (prt_level.ge.1) then
1963	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1964	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1965	endif
1966	c
1967	c Appeler la convection (au choix)
1968	c
1969	DO k = 1, klev
1970	DO i = 1, klon
1971	conv_q(i,k) = d_q_dyn(i,k)
1972	. + d_q_vdf(i,k)/dtime
1973	conv_t(i,k) = d_t_dyn(i,k)
1974	. + d_t_vdf(i,k)/dtime
1975	ENDDO
1976	ENDDO
1977	IF (check) THEN
1978	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1979	WRITE(lunout,*) "avantcon=", za
1980	ENDIF
1981	zx_ajustq = .FALSE.
1982	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1983	IF (zx_ajustq) THEN
1984	DO i = 1, klon
1985	z_avant(i) = 0.0
1986	ENDDO
1987	DO k = 1, klev
1988	DO i = 1, klon
1989	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1990	. *(paprs(i,k)-paprs(i,k+1))/RG
1991	ENDDO
1992	ENDDO
1993	ENDIF
1994
1995	c Calcule de vitesse verticale a partir de flux de masse verticale
1996	DO k = 1, klev
1997	DO i = 1, klon
1998	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1999	END DO
2000	END DO
2001	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ',
2002	$ omega(igout, :)
2003
2004	IF (iflag_con.EQ.1) THEN
2005	stop'reactiver le call conlmd dans physiq.F'
2006	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2007	c . d_t_con, d_q_con,
2008	c . rain_con, snow_con, ibas_con, itop_con)
2009	ELSE IF (iflag_con.EQ.2) THEN
2010	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2011	e conv_t, conv_q, -evap, omega,
2012	s d_t_con, d_q_con, rain_con, snow_con,
2013	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2014	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2015	d_u_con = 0.
2016	d_v_con = 0.
2017
2018	WHERE (rain_con < 0.) rain_con = 0.
2019	WHERE (snow_con < 0.) snow_con = 0.
2020	DO i = 1, klon
2021	ibas_con(i) = klev+1 - kcbot(i)
2022	itop_con(i) = klev+1 - kctop(i)
2023	ENDDO
2024	ELSE IF (iflag_con.GE.3) THEN
2025	c nb of tracers for the KE convection:
2026	c MAF la partie traceurs est faite dans phytrac
2027	c on met ntra=1 pour limiter les appels mais on peut
2028	c supprimer les calculs / ftra.
2029	ntra = 1
2030
2031	c=====================================================================================
2032	cajout pour la parametrisation des poches froides:
2033	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2034	do k=1,klev
2035	do i=1,klon
2036	if (iflag_wake.eq.1) then
2037	t_wake(i,k) = t_seri(i,k)
2038	. +(1-wake_s(i))*wake_deltat(i,k)
2039	q_wake(i,k) = q_seri(i,k)
2040	. +(1-wake_s(i))*wake_deltaq(i,k)
2041	t_undi(i,k) = t_seri(i,k)
2042	. -wake_s(i)*wake_deltat(i,k)
2043	q_undi(i,k) = q_seri(i,k)
2044	. -wake_s(i)*wake_deltaq(i,k)
2045	else
2046	t_wake(i,k) = t_seri(i,k)
2047	q_wake(i,k) = q_seri(i,k)
2048	t_undi(i,k) = t_seri(i,k)
2049	q_undi(i,k) = q_seri(i,k)
2050	endif
2051	enddo
2052	enddo
2053
2054	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2055	cc-- pour le soulevement des particules dans le modele convectif
2056	c
2057	do i = 1,klon
2058	ALE(i) = 0.
2059	ALP(i) = 0.
2060	enddo
2061	c
2062	ccalcul de ale_wake et alp_wake
2063	do i = 1,klon
2064	if (iflag_wake.eq.1) then
2065	ale_wake(i) = 0.5wake_cstar(i)*2
2066	alp_wake(i) = wake_fip(i)
2067	else
2068	ale_wake(i) = 0.
2069	alp_wake(i) = 0.
2070	endif
2071	enddo
2072	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2073	cdans le thermique sinon
2074	if (iflag_coupl.eq.0) then
2075	if (debut) print*,'ALE et ALP imposes'
2076	do i = 1,klon
2077	con ne couple que ale
2078	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2079	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2080	con ne couple que alp
2081	c ALP(i) = alp_wake(i) + Alp_bl(i)
2082	ALP(i) = alp_wake(i) + alp_bl_prescr
2083	enddo
2084	else
2085	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2086	do i = 1,klon
2087	ALE(i) = max(ale_wake(i),Ale_bl(i))
2088	ALP(i) = alp_wake(i) + Alp_bl(i)
2089	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2090	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2091	enddo
2092	endif
2093	do i=1,klon
2094	if (alp(i)>alp_max) then
2095	IF(prt_level>9)WRITE(lunout,*) &
2096	& 'WARNING SUPER ALP (seuil=',alp_max,
2097	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2098	alp(i)=alp_max
2099	endif
2100	if (ale(i)>ale_max) then
2101	IF(prt_level>9)WRITE(lunout,*) &
2102	& 'WARNING SUPER ALE (seuil=',ale_max,
2103	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2104	ale(i)=ale_max
2105	endif
2106	enddo
2107
2108	cfin calcul ale et alp
2109	c=================================================================================================
2110
2111
2112	c sb, oct02:
2113	c Schema de convection modularise et vectorise:
2114	c (driver commun aux versions 3 et 4)
2115	c
2116	IF (ok_cvl) THEN ! new driver for convectL
2117
2118	CALL concvl (iflag_con,iflag_clos,
2119	. dtime,paprs,pplay,t_undi,q_undi,
2120	. t_wake,q_wake,wake_s,
2121	. u_seri,v_seri,tr_seri,nbtr,
2122	. ALE,ALP,
2123	. ema_work1,ema_work2,
2124	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2125	. rain_con, snow_con, ibas_con, itop_con, sigd,
2126	. upwd,dnwd,dnwd0,
2127	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2128	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2129	. pmflxr,pmflxs,da,phi,mp,
2130	. ftd,fqd,lalim_conv,wght_th)
2131
2132	cIM begin
2133	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2134	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2135	cIM end
2136	cIM cf. FH
2137	clwcon0=qcondc
2138	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2139
2140	do i = 1, klon
2141	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2142	enddo
2143
2144	ELSE ! ok_cvl
2145	c MAF conema3 ne contient pas les traceurs
2146	CALL conema3 (dtime,
2147	. paprs,pplay,t_seri,q_seri,
2148	. u_seri,v_seri,tr_seri,ntra,
2149	. ema_work1,ema_work2,
2150	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2151	. rain_con, snow_con, ibas_con, itop_con,
2152	. upwd,dnwd,dnwd0,bas,top,
2153	. Ma,cape,tvp,rflag,
2154	. pbase
2155	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2156	. ,clwcon0)
2157
2158	ENDIF ! ok_cvl
2159
2160	c
2161	c Correction precip
2162	rain_con = rain_con * cvl_corr
2163	snow_con = snow_con * cvl_corr
2164	c
2165
2166	IF (.NOT. ok_gust) THEN
2167	do i = 1, klon
2168	wd(i)=0.0
2169	enddo
2170	ENDIF
2171
2172	c =================================================================== c
2173	c Calcul des proprietes des nuages convectifs
2174	c
2175
2176	c calcul des proprietes des nuages convectifs
2177	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2178	call clouds_gno
2179	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2180
2181	c =================================================================== c
2182
2183	DO i = 1, klon
2184	ema_pcb(i) = paprs(i,ibas_con(i))
2185	ENDDO
2186	DO i = 1, klon
2187	! L'idicage de itop_con peut cacher un pb potentiel
2188	! FH sous la dictee de JYG, CR
2189	ema_pct(i) = paprs(i,itop_con(i)+1)
2190
2191	if (itop_con(i).gt.klev-3) then
2192	if(prt_level >= 9) then
2193	write(lunout,*)'La convection monte trop haut '
2194	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2195	endif
2196	endif
2197	ENDDO
2198	DO i = 1, klon
2199	ema_cbmf(i) = ema_workcbmf(i)
2200	ENDDO
2201	ELSE IF (iflag_con.eq.0) THEN
2202	write(lunout,*) 'On n appelle pas la convection'
2203	clwcon0=0.
2204	rnebcon0=0.
2205	d_t_con=0.
2206	d_q_con=0.
2207	d_u_con=0.
2208	d_v_con=0.
2209	rain_con=0.
2210	snow_con=0.
2211	bas=1
2212	top=1
2213	ELSE
2214	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2215	CALL abort
2216	ENDIF
2217
2218	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2219	c . d_u_con, d_v_con)
2220
2221	!-----------------------------------------------------------------------------------------
2222	! ajout des tendances de la diffusion turbulente
2223	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2224	!-----------------------------------------------------------------------------------------
2225
2226	if (mydebug) then
2227	call writefield_phy('u_seri',u_seri,llm)
2228	call writefield_phy('v_seri',v_seri,llm)
2229	call writefield_phy('t_seri',t_seri,llm)
2230	call writefield_phy('q_seri',q_seri,llm)
2231	endif
2232
2233	cIM
2234	IF (ip_ebil_phy.ge.2) THEN
2235	ztit='after convect'
2236	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2237	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2238	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2239	call diagphy(airephy,ztit,ip_ebil_phy
2240	e , zero_v, zero_v, zero_v, zero_v, zero_v
2241	e , zero_v, rain_con, snow_con, ztsol
2242	e , d_h_vcol, d_qt, d_ec
2243	s , fs_bound, fq_bound )
2244	END IF
2245	C
2246	IF (check) THEN
2247	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2248	WRITE(lunout,*)"aprescon=", za
2249	zx_t = 0.0
2250	za = 0.0
2251	DO i = 1, klon
2252	za = za + airephy(i)/FLOAT(klon)
2253	zx_t = zx_t + (rain_con(i)+
2254	. snow_con(i))*airephy(i)/FLOAT(klon)
2255	ENDDO
2256	zx_t = zx_t/za*dtime
2257	WRITE(lunout,*)"Precip=", zx_t
2258	ENDIF
2259	IF (zx_ajustq) THEN
2260	DO i = 1, klon
2261	z_apres(i) = 0.0
2262	ENDDO
2263	DO k = 1, klev
2264	DO i = 1, klon
2265	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2266	. *(paprs(i,k)-paprs(i,k+1))/RG
2267	ENDDO
2268	ENDDO
2269	DO i = 1, klon
2270	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2271	. /z_apres(i)
2272	ENDDO
2273	DO k = 1, klev
2274	DO i = 1, klon
2275	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2276	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2277	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2278	ENDIF
2279	ENDDO
2280	ENDDO
2281	ENDIF
2282	zx_ajustq=.FALSE.
2283
2284	c
2285	c=============================================================================
2286	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2287	cpour la couche limite diffuse pour l instant
2288	c
2289	if (iflag_wake.eq.1) then
2290	DO k=1,klev
2291	DO i=1,klon
2292	dt_dwn(i,k) = ftd(i,k)
2293	wdt_PBL(i,k) = 0.
2294	dq_dwn(i,k) = fqd(i,k)
2295	wdq_PBL(i,k) = 0.
2296	M_dwn(i,k) = dnwd0(i,k)
2297	M_up(i,k) = upwd(i,k)
2298	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2299	udt_PBL(i,k) = 0.
2300	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2301	udq_PBL(i,k) = 0.
2302	ENDDO
2303	ENDDO
2304	c
2305	ccalcul caracteristiques de la poche froide
2306	call calWAKE (paprs,pplay,dtime
2307	: ,t_seri,q_seri,omega
2308	: ,dt_dwn,dq_dwn,M_dwn,M_up
2309	: ,dt_a,dq_a,sigd
2310	: ,wdt_PBL,wdq_PBL
2311	: ,udt_PBL,udq_PBL
2312	o ,wake_deltat,wake_deltaq,wake_dth
2313	o ,wake_h,wake_s,wake_dens
2314	o ,wake_pe,wake_fip,wake_gfl
2315	o ,dt_wake,dq_wake
2316	o ,wake_k, t_undi,q_undi
2317	o ,wake_omgbdth,wake_dp_omgb
2318	o ,wake_dtKE,wake_dqKE
2319	o ,wake_dtPBL,wake_dqPBL
2320	o ,wake_omg,wake_dp_deltomg
2321	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2322	o ,wake_ddeltat,wake_ddeltaq)
2323	c
2324	!-----------------------------------------------------------------------------------------
2325	! ajout des tendances des poches froides
2326	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2327	! coherent avec les autres d_t_...
2328	d_t_wake(:,:)=dt_wake(:,:)*dtime
2329	d_q_wake(:,:)=dq_wake(:,:)*dtime
2330	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2331	!-----------------------------------------------------------------------------------------
2332
2333	endif
2334	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2335	c
2336	c===================================================================
2337	c Convection seche (thermiques ou ajustement)
2338	c===================================================================
2339	c
2340	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2341	s ,seuil_inversion,weak_inversion,dthmin)
2342
2343
2344
2345	d_t_ajsb(:,:)=0.
2346	d_q_ajsb(:,:)=0.
2347	d_t_ajs(:,:)=0.
2348	d_u_ajs(:,:)=0.
2349	d_v_ajs(:,:)=0.
2350	d_q_ajs(:,:)=0.
2351	clwcon0th(:,:)=0.
2352	c
2353	fm_therm(:,:)=0.
2354	entr_therm(:,:)=0.
2355	detr_therm(:,:)=0.
2356	c
2357	IF(prt_level>9)WRITE(lunout,*)
2358	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2359	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2360	if(iflag_thermals.lt.0) then
2361	c Rien
2362	c ====
2363	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2364
2365
2366	else
2367
2368	c Thermiques
2369	c ==========
2370	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2371	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2372
2373
2374	if (iflag_thermals.gt.1) then
2375	call calltherm(pdtphys
2376	s ,pplay,paprs,pphi,weak_inversion
2377	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2378	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2379	s ,fm_therm,entr_therm,detr_therm
2380	s ,zqasc,clwcon0th,lmax_th,ratqscth
2381	s ,ratqsdiff,zqsatth
2382	con rajoute ale et alp, et les caracteristiques de la couche alim
2383	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
2384	endif
2385
2386
2387	c Ajustement sec
2388	c ==============
2389
2390	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2391	! a partir du sommet des thermiques.
2392	! Dans le cas contraire, on demarre au niveau 1.
2393
2394	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2395
2396	if(iflag_thermals.eq.0) then
2397	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2398	limbas(:)=1
2399	else
2400	limbas(:)=lmax_th(:)
2401	endif
2402
2403	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2404	! pour des test de convergence numerique.
2405	! Le nouveau ajsec est a priori mieux, meme pour le cas
2406	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2407	! non nulles numeriquement pour des mailles non concernees.
2408
2409	if (iflag_thermals.eq.0) then
2410	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2411	s , d_t_ajsb, d_q_ajsb)
2412	else
2413	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2414	s , d_t_ajsb, d_q_ajsb)
2415	endif
2416
2417	!-----------------------------------------------------------------------------------------
2418	! ajout des tendances de l'ajustement sec ou des thermiques
2419	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2420	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2421	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2422
2423	!-----------------------------------------------------------------------------------------
2424
2425	endif
2426
2427	endif
2428	c
2429	c===================================================================
2430	cIM
2431	IF (ip_ebil_phy.ge.2) THEN
2432	ztit='after dry_adjust'
2433	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2434	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2435	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2436	END IF
2437
2438
2439	c-------------------------------------------------------------------------
2440	c Caclul des ratqs
2441	c-------------------------------------------------------------------------
2442
2443	c print*,'calcul des ratqs'
2444	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2445	c ----------------
2446	c on ecrase le tableau ratqsc calcule par clouds_gno
2447	if (iflag_cldcon.eq.1) then
2448	do k=1,klev
2449	do i=1,klon
2450	if(ptconv(i,k)) then
2451	ratqsc(i,k)=ratqsbas
2452	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2453	else
2454	ratqsc(i,k)=0.
2455	endif
2456	enddo
2457	enddo
2458
2459	c-----------------------------------------------------------------------
2460	c par nversion de la fonction log normale
2461	c-----------------------------------------------------------------------
2462	else if (iflag_cldcon.eq.4) then
2463	ptconvth(:,:)=.false.
2464	ratqsc(:,:)=0.
2465	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2466	call clouds_gno
2467	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2468	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2469
2470	c-----------------------------------------------------------------------
2471	c par calcul direct de l'ecart-type
2472	c-----------------------------------------------------------------------
2473
2474	else if (iflag_cldcon>=5) then
2475	wmax_th(:)=0.
2476	zmax_th(:)=0.
2477	do k=1,klev
2478	do i=1,klon
2479	wmax_th(i)=max(wmax_th(i),zw2(i,k))
2480	if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
2481	enddo
2482	enddo
2483	tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
2484	print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
2485
2486	c On impose que l'air autour de la fraction couverte par le thermique
2487	c plus son air detraine durant tau_overturning_th soit superieur
2488	c a 0.1 q_seri
2489	zz=0.1
2490	do k=1,klev
2491	do i=1,klon
2492	lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
2493	s tau_overturning_th(i)*detr_therm(i,k)
2494	s *rg/(paprs(i,k)-paprs(i,k+1))
2495	znum=(1.-zz)*q_seri(i,k)
2496	zden=zqasc(i,k)-zz*q_seri(i,k)
2497	if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
2498	lambda_th(i,k)=min(lambda_th(i,k),0.9)
2499	enddo
2500	enddo
2501
2502	if(iflag_cldcon==5) then
2503	do k=1,klev
2504	do i=1,klon
2505	ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
2506	s abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
2507	enddo
2508	enddo
2509	else if(iflag_cldcon==6) then
2510	do k=1,klev
2511	do i=1,klon
2512	ratqsc(i,k)=sqrt(lambda_th(i,k))*
2513	s (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
2514	enddo
2515	enddo
2516	endif
2517
2518	endif
2519
2520	c ratqs stables
2521	c -------------
2522
2523	if (iflag_ratqs.eq.0) then
2524
2525	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2526	do k=1,klev
2527	do i=1, klon
2528	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2529	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2530	enddo
2531	enddo
2532
2533	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2534	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2535	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2536	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2537	! Il s'agit de differents tests dans la phase de reglage du modele
2538	! avec thermiques.
2539
2540	else if (iflag_ratqs.eq.1) then
2541
2542	do k=1,klev
2543	do i=1, klon
2544	if (pplay(i,k).ge.60000.) then
2545	ratqss(i,k)=ratqsbas
2546	else if ((pplay(i,k).ge.30000.).and.
2547	s (pplay(i,k).lt.60000.)) then
2548	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2549	s (60000.-pplay(i,k))/(60000.-30000.)
2550	else
2551	ratqss(i,k)=ratqshaut
2552	endif
2553	enddo
2554	enddo
2555
2556	else
2557
2558	do k=1,klev
2559	do i=1, klon
2560	if (pplay(i,k).ge.60000.) then
2561	ratqss(i,k)=ratqsbas
2562	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2563	else if ((pplay(i,k).ge.30000.).and.
2564	s (pplay(i,k).lt.60000.)) then
2565	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2566	s (60000.-pplay(i,k))/(60000.-30000.)
2567	else
2568	ratqss(i,k)=ratqshaut
2569	endif
2570	enddo
2571	enddo
2572	endif
2573
2574
2575
2576
2577	c ratqs final
2578	c -----------
2579
2580	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2581	s .or.iflag_cldcon.ge.4) then
2582
2583	! On ajoute une constante au ratqsc*2 pour tenir compte de
2584	! fluctuations turbulentes de petite echelle
2585
2586	do k=1,klev
2587	do i=1,klon
2588	if ((fm_therm(i,k).gt.1.e-10)) then
2589	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2590	endif
2591	enddo
2592	enddo
2593
2594	! les ratqs sont une combinaison de ratqss et ratqsc
2595	! print*,'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
2596
2597	if (tau_ratqs>1.e-10) then
2598	facteur=exp(-pdtphys/tau_ratqs)
2599	else
2600	facteur=0.
2601	endif
2602	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2603	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2604	! FH 22/09/2009
2605	! La ligne ci-dessous faisait osciller le modele et donnait une solution
2606	! assymptotique bidon et dÃ©pendant fortement du pas de temps.
2607	! ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2608	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2609	ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
2610	else
2611	! on ne prend que le ratqs stable pour fisrtilp
2612	ratqs(:,:)=ratqss(:,:)
2613	endif
2614
2615
2616	c
2617	c Appeler le processus de condensation a grande echelle
2618	c et le processus de precipitation
2619	c-------------------------------------------------------------------------
2620	CALL fisrtilp(dtime,paprs,pplay,
2621	. t_seri, q_seri,ptconv,ratqs,
2622	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2623	. rain_lsc, snow_lsc,
2624	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2625	. frac_impa, frac_nucl,
2626	. prfl, psfl, rhcl)
2627
2628	WHERE (rain_lsc < 0) rain_lsc = 0.
2629	WHERE (snow_lsc < 0) snow_lsc = 0.
2630	!-----------------------------------------------------------------------------------------
2631	! ajout des tendances de la diffusion turbulente
2632	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2633	!-----------------------------------------------------------------------------------------
2634	DO k = 1, klev
2635	DO i = 1, klon
2636	cldfra(i,k) = rneb(i,k)
2637	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2638	ENDDO
2639	ENDDO
2640	IF (check) THEN
2641	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2642	WRITE(lunout,*)"apresilp=", za
2643	zx_t = 0.0
2644	za = 0.0
2645	DO i = 1, klon
2646	za = za + airephy(i)/FLOAT(klon)
2647	zx_t = zx_t + (rain_lsc(i)
2648	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2649	ENDDO
2650	zx_t = zx_t/za*dtime
2651	WRITE(lunout,*)"Precip=", zx_t
2652	ENDIF
2653	cIM
2654	IF (ip_ebil_phy.ge.2) THEN
2655	ztit='after fisrt'
2656	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2657	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2658	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2659	call diagphy(airephy,ztit,ip_ebil_phy
2660	e , zero_v, zero_v, zero_v, zero_v, zero_v
2661	e , zero_v, rain_lsc, snow_lsc, ztsol
2662	e , d_h_vcol, d_qt, d_ec
2663	s , fs_bound, fq_bound )
2664	END IF
2665
2666	if (mydebug) then
2667	call writefield_phy('u_seri',u_seri,llm)
2668	call writefield_phy('v_seri',v_seri,llm)
2669	call writefield_phy('t_seri',t_seri,llm)
2670	call writefield_phy('q_seri',q_seri,llm)
2671	endif
2672
2673	c
2674	c-------------------------------------------------------------------
2675	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2676	c-------------------------------------------------------------------
2677
2678	c 1. NUAGES CONVECTIFS
2679	c
2680	cIM cf FH
2681	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2682	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2683	snow_tiedtke=0.
2684	c print*,'avant calcul de la pseudo precip '
2685	c print*,'iflag_cldcon',iflag_cldcon
2686	if (iflag_cldcon.eq.-1) then
2687	rain_tiedtke=rain_con
2688	else
2689	c print*,'calcul de la pseudo precip '
2690	rain_tiedtke=0.
2691	c print*,'calcul de la pseudo precip 0'
2692	do k=1,klev
2693	do i=1,klon
2694	if (d_q_con(i,k).lt.0.) then
2695	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2696	s *(paprs(i,k)-paprs(i,k+1))/rg
2697	endif
2698	enddo
2699	enddo
2700	endif
2701	c
2702	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2703	c
2704
2705	c Nuages diagnostiques pour Tiedtke
2706	CALL diagcld1(paprs,pplay,
2707	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2708	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2709	. diafra,dialiq)
2710	DO k = 1, klev
2711	DO i = 1, klon
2712	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2713	cldliq(i,k) = dialiq(i,k)
2714	cldfra(i,k) = diafra(i,k)
2715	ENDIF
2716	ENDDO
2717	ENDDO
2718
2719	ELSE IF (iflag_cldcon.ge.3) THEN
2720	c On prend pour les nuages convectifs le max du calcul de la
2721	c convection et du calcul du pas de temps precedent diminue d'un facteur
2722	c facttemps
2723	facteur = pdtphys *facttemps
2724	do k=1,klev
2725	do i=1,klon
2726	rnebcon(i,k)=rnebcon(i,k)*facteur
2727	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2728	s then
2729	rnebcon(i,k)=rnebcon0(i,k)
2730	clwcon(i,k)=clwcon0(i,k)
2731	endif
2732	enddo
2733	enddo
2734
2735	c
2736	cjq - introduce the aerosol direct and first indirect radiative forcings
2737	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2738	IF (ok_ade.OR.ok_aie) THEN
2739	IF (.NOT. aerosol_couple)
2740	& CALL readaerosol_optic(
2741	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref,
2742	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs,
2743	& mass_solu_aero, mass_solu_aero_pi,
2744	& tau_aero, piz_aero, cg_aero,
2745	& tausum_aero, tau3d_aero)
2746	ELSE
2747	tau_aero(:,:,:,:) = 0.
2748	piz_aero(:,:,:,:) = 0.
2749	cg_aero(:,:,:,:) = 0.
2750	ENDIF
2751
2752	cIM calcul nuages par le simulateur ISCCP
2753	c
2754	#ifdef histISCCP
2755	IF (ok_isccp) THEN
2756	c
2757	cIM lecture invtau, tautab des fichiers formattes
2758	c
2759	IF (debut) THEN
2760	c$OMP MASTER
2761	c
2762	open(99,file='tautab.formatted', FORM='FORMATTED')
2763	read(99,'(f30.20)') tautab_omp
2764	close(99)
2765	c
2766	open(99,file='invtau.formatted',form='FORMATTED')
2767	read(99,'(i10)') invtau_omp
2768
2769	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
2770	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2771
2772	close(99)
2773	c$OMP END MASTER
2774	c$OMP BARRIER
2775	tautab=tautab_omp
2776	invtau=invtau_omp
2777	c
2778	ENDIF !debut
2779	c
2780	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2781	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2782	#include "calcul_simulISCCP.h"
2783	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2784	ENDIF !ok_isccp
2785	#endif
2786
2787	c On prend la somme des fractions nuageuses et des contenus en eau
2788	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2789	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2790
2791	ENDIF
2792	c
2793	c 2. NUAGES STARTIFORMES
2794	c
2795	IF (ok_stratus) THEN
2796	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2797	DO k = 1, klev
2798	DO i = 1, klon
2799	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2800	cldliq(i,k) = dialiq(i,k)
2801	cldfra(i,k) = diafra(i,k)
2802	ENDIF
2803	ENDDO
2804	ENDDO
2805	ENDIF
2806	c
2807	c Precipitation totale
2808	c
2809	DO i = 1, klon
2810	rain_fall(i) = rain_con(i) + rain_lsc(i)
2811	snow_fall(i) = snow_con(i) + snow_lsc(i)
2812	ENDDO
2813	cIM
2814	IF (ip_ebil_phy.ge.2) THEN
2815	ztit="after diagcld"
2816	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2817	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2818	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2819	END IF
2820	c
2821	c Calculer l'humidite relative pour diagnostique
2822	c
2823	DO k = 1, klev
2824	DO i = 1, klon
2825	zx_t = t_seri(i,k)
2826	IF (thermcep) THEN
2827	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2828	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2829	zx_qs = MIN(0.5,zx_qs)
2830	zcor = 1./(1.-retv*zx_qs)
2831	zx_qs = zx_qs*zcor
2832	ELSE
2833	IF (zx_t.LT.t_coup) THEN
2834	zx_qs = qsats(zx_t)/pplay(i,k)
2835	ELSE
2836	zx_qs = qsatl(zx_t)/pplay(i,k)
2837	ENDIF
2838	ENDIF
2839	zx_rh(i,k) = q_seri(i,k)/zx_qs
2840	zqsat(i,k)=zx_qs
2841	ENDDO
2842	ENDDO
2843
2844	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2845	c equivalente a 2m (tpote) pour diagnostique
2846	c
2847	DO i = 1, klon
2848	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2849	IF (thermcep) THEN
2850	IF(zt2m(i).LT.RTT) then
2851	Lheat=RLSTT
2852	ELSE
2853	Lheat=RLVTT
2854	ENDIF
2855	ELSE
2856	IF (zt2m(i).LT.RTT) THEN
2857	Lheat=RLSTT
2858	ELSE
2859	Lheat=RLVTT
2860	ENDIF
2861	ENDIF
2862	tpote(i) = tpot(i)*
2863	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2864	ENDDO
2865
2866	IF (config_inca /= 'none') THEN
2867	#ifdef INCA
2868	CALL VTe(VTphysiq)
2869	CALL VTb(VTinca)
2870	calday = FLOAT(days_elapsed + 1) + jH_cur
2871
2872	call chemtime(itap+itau_phy-1, date0, dtime)
2873	IF (config_inca == 'aero') THEN
2874	CALL AEROSOL_METEO_CALC(
2875	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2876	$ prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
2877	END IF
2878
2879	zxsnow_dummy(:) = 0.0
2880
2881	CALL chemhook_begin (calday,
2882	$ days_elapsed+1,
2883	$ jH_cur,
2884	$ pctsrf(1,1),
2885	$ rlat,
2886	$ rlon,
2887	$ airephy,
2888	$ paprs,
2889	$ pplay,
2890	$ coefh,
2891	$ pphi,
2892	$ t_seri,
2893	$ u,
2894	$ v,
2895	$ wo(:, :, 1),
2896	$ q_seri,
2897	$ zxtsol,
2898	$ zxsnow_dummy,
2899	$ solsw,
2900	$ albsol1,
2901	$ rain_fall,
2902	$ snow_fall,
2903	$ itop_con,
2904	$ ibas_con,
2905	$ cldfra,
2906	$ iim,
2907	$ jjm,
2908	$ tr_seri,
2909	$ ftsol,
2910	$ paprs,
2911	$ cdragh,
2912	$ cdragm,
2913	$ pctsrf,
2914	$ pdtphys,
2915	$ itap)
2916
2917	CALL VTe(VTinca)
2918	CALL VTb(VTphysiq)
2919	#endif
2920	END IF !config_inca /= 'none'
2921	c
2922	c Calculer les parametres optiques des nuages et quelques
2923	c parametres pour diagnostiques:
2924	c
2925
2926	IF (aerosol_couple) THEN
2927	mass_solu_aero(:,:) = ccm(:,:,1)
2928	mass_solu_aero_pi(:,:) = ccm(:,:,2)
2929	END IF
2930
2931	if (ok_newmicro) then
2932	CALL newmicro (paprs, pplay,ok_newmicro,
2933	. t_seri, cldliq, cldfra, cldtau, cldemi,
2934	. cldh, cldl, cldm, cldt, cldq,
2935	. flwp, fiwp, flwc, fiwc,
2936	e ok_aie,
2937	e mass_solu_aero, mass_solu_aero_pi,
2938	e bl95_b0, bl95_b1,
2939	s cldtaupi, re, fl, ref_liq, ref_ice)
2940	else
2941	CALL nuage (paprs, pplay,
2942	. t_seri, cldliq, cldfra, cldtau, cldemi,
2943	. cldh, cldl, cldm, cldt, cldq,
2944	e ok_aie,
2945	e mass_solu_aero, mass_solu_aero_pi,
2946	e bl95_b0, bl95_b1,
2947	s cldtaupi, re, fl)
2948
2949	endif
2950	c
2951	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2952	c
2953	IF (MOD(itaprad,radpas).EQ.0) THEN
2954
2955	DO i = 1, klon
2956	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2957	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2958	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2959	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2960	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2961	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2962	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2963	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2964	ENDDO
2965
2966	if (mydebug) then
2967	call writefield_phy('u_seri',u_seri,llm)
2968	call writefield_phy('v_seri',v_seri,llm)
2969	call writefield_phy('t_seri',t_seri,llm)
2970	call writefield_phy('q_seri',q_seri,llm)
2971	endif
2972
2973	IF (aerosol_couple) THEN
2974	#ifdef INCA
2975	CALL radlwsw_inca
2976	e (kdlon,kflev,dist, rmu0, fract, solaire,
2977	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2978	e wo(:, :, 1),
2979	e cldfra, cldemi, cldtau,
2980	s heat,heat0,cool,cool0,radsol,albpla,
2981	s topsw,toplw,solsw,sollw,
2982	s sollwdown,
2983	s topsw0,toplw0,solsw0,sollw0,
2984	s lwdn0, lwdn, lwup0, lwup,
2985	s swdn0, swdn, swup0, swup,
2986	e ok_ade, ok_aie,
2987	e tau_aero, piz_aero, cg_aero,
2988	s topswad_aero, solswad_aero,
2989	s topswad0_aero, solswad0_aero,
2990	s topsw_aero, topsw0_aero,
2991	s solsw_aero, solsw0_aero,
2992	e cldtaupi,
2993	s topswai_aero, solswai_aero)
2994
2995	#endif
2996	ELSE
2997
2998	CALL radlwsw
2999	e (dist, rmu0, fract,
3000	e paprs, pplay,zxtsol,albsol1, albsol2,
3001	e t_seri,q_seri,wo,
3002	e cldfra, cldemi, cldtau,
3003	e ok_ade, ok_aie,
3004	e tau_aero, piz_aero, cg_aero,
3005	e cldtaupi,new_aod,
3006	e zqsat, flwc, fiwc,
3007	s heat,heat0,cool,cool0,radsol,albpla,
3008	s topsw,toplw,solsw,sollw,
3009	s sollwdown,
3010	s topsw0,toplw0,solsw0,sollw0,
3011	s lwdn0, lwdn, lwup0, lwup,
3012	s swdn0, swdn, swup0, swup,
3013	s topswad_aero, solswad_aero,
3014	s topswai_aero, solswai_aero,
3015	o topswad0_aero, solswad0_aero,
3016	o topsw_aero, topsw0_aero,
3017	o solsw_aero, solsw0_aero,
3018	o topswcf_aero, solswcf_aero)
3019
3020
3021	ENDIF ! aerosol_couple
3022	itaprad = 0
3023	ENDIF ! MOD(itaprad,radpas)
3024	itaprad = itaprad + 1
3025
3026	if (iflag_radia.eq.0) then
3027	print *,'--------------------------------------------------'
3028	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3029	print *,'>>>> heat et cool mis a zero '
3030	print *,'--------------------------------------------------'
3031	heat=0.
3032	cool=0.
3033	endif
3034
3035	c
3036	c Ajouter la tendance des rayonnements (tous les pas)
3037	c
3038	DO k = 1, klev
3039	DO i = 1, klon
3040	t_seri(i,k) = t_seri(i,k)
3041	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
3042	ENDDO
3043	ENDDO
3044	c
3045	if (mydebug) then
3046	call writefield_phy('u_seri',u_seri,llm)
3047	call writefield_phy('v_seri',v_seri,llm)
3048	call writefield_phy('t_seri',t_seri,llm)
3049	call writefield_phy('q_seri',q_seri,llm)
3050	endif
3051
3052	cIM
3053	IF (ip_ebil_phy.ge.2) THEN
3054	ztit='after rad'
3055	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3056	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3057	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3058	call diagphy(airephy,ztit,ip_ebil_phy
3059	e , topsw, toplw, solsw, sollw, zero_v
3060	e , zero_v, zero_v, zero_v, ztsol
3061	e , d_h_vcol, d_qt, d_ec
3062	s , fs_bound, fq_bound )
3063	END IF
3064	c
3065	c
3066	c Calculer l'hydrologie de la surface
3067	c
3068	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3069	c . agesno, ftsol,fqsurf,fsnow, ruis)
3070	c
3071
3072	c
3073	c Calculer le bilan du sol et la derive de temperature (couplage)
3074	c
3075	DO i = 1, klon
3076	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3077	c a la demande de JLD
3078	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3079	ENDDO
3080	c
3081	cmoddeblott(jan95)
3082	c Appeler le programme de parametrisation de l'orographie
3083	c a l'echelle sous-maille:
3084	c
3085	IF (ok_orodr) THEN
3086	c
3087	c selection des points pour lesquels le shema est actif:
3088	igwd=0
3089	DO i=1,klon
3090	itest(i)=0
3091	c IF ((zstd(i).gt.10.0)) THEN
3092	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3093	itest(i)=1
3094	igwd=igwd+1
3095	idx(igwd)=i
3096	ENDIF
3097	ENDDO
3098	c igwdim=MAX(1,igwd)
3099	c
3100	IF (ok_strato) THEN
3101
3102	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3103	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3104	e igwd,idx,itest,
3105	e t_seri, u_seri, v_seri,
3106	s zulow, zvlow, zustrdr, zvstrdr,
3107	s d_t_oro, d_u_oro, d_v_oro)
3108
3109	ELSE
3110	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3111	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3112	e igwd,idx,itest,
3113	e t_seri, u_seri, v_seri,
3114	s zulow, zvlow, zustrdr, zvstrdr,
3115	s d_t_oro, d_u_oro, d_v_oro)
3116	ENDIF
3117	c
3118	c ajout des tendances
3119	!-----------------------------------------------------------------------------------------
3120	! ajout des tendances de la trainee de l'orographie
3121	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3122	!-----------------------------------------------------------------------------------------
3123	c
3124	ENDIF ! fin de test sur ok_orodr
3125	c
3126	if (mydebug) then
3127	call writefield_phy('u_seri',u_seri,llm)
3128	call writefield_phy('v_seri',v_seri,llm)
3129	call writefield_phy('t_seri',t_seri,llm)
3130	call writefield_phy('q_seri',q_seri,llm)
3131	endif
3132
3133	IF (ok_orolf) THEN
3134	c
3135	c selection des points pour lesquels le shema est actif:
3136	igwd=0
3137	DO i=1,klon
3138	itest(i)=0
3139	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3140	itest(i)=1
3141	igwd=igwd+1
3142	idx(igwd)=i
3143	ENDIF
3144	ENDDO
3145	c igwdim=MAX(1,igwd)
3146	c
3147	IF (ok_strato) THEN
3148
3149	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3150	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3151	e igwd,idx,itest,
3152	e t_seri, u_seri, v_seri,
3153	s zulow, zvlow, zustrli, zvstrli,
3154	s d_t_lif, d_u_lif, d_v_lif )
3155
3156	ELSE
3157	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3158	e rlat,zmea,zstd,zpic,
3159	e itest,
3160	e t_seri, u_seri, v_seri,
3161	s zulow, zvlow, zustrli, zvstrli,
3162	s d_t_lif, d_u_lif, d_v_lif)
3163	ENDIF
3164	c
3165	!-----------------------------------------------------------------------------------------
3166	! ajout des tendances de la portance de l'orographie
3167	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3168	!-----------------------------------------------------------------------------------------
3169	c
3170	ENDIF ! fin de test sur ok_orolf
3171	C HINES GWD PARAMETRIZATION
3172
3173	IF (ok_hines) then
3174
3175	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3176	i rlat,t_seri,u_seri,v_seri,
3177	o zustrhi,zvstrhi,
3178	o d_t_hin, d_u_hin, d_v_hin)
3179	c
3180	c ajout des tendances
3181	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3182
3183	ENDIF
3184	c
3185
3186	c
3187	cIM cf. FLott BEG
3188	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3189
3190	if (mydebug) then
3191	call writefield_phy('u_seri',u_seri,llm)
3192	call writefield_phy('v_seri',v_seri,llm)
3193	call writefield_phy('t_seri',t_seri,llm)
3194	call writefield_phy('q_seri',q_seri,llm)
3195	endif
3196
3197	DO i = 1, klon
3198	zustrph(i)=0.
3199	zvstrph(i)=0.
3200	ENDDO
3201	DO k = 1, klev
3202	DO i = 1, klon
3203	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3204	c (paprs(i,k)-paprs(i,k+1))/rg
3205	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3206	c (paprs(i,k)-paprs(i,k+1))/rg
3207	ENDDO
3208	ENDDO
3209	c
3210	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3211	c
3212	IF (is_sequential) THEN
3213
3214	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
3215	C ra,rg,romega,
3216	C rlat,rlon,pphis,
3217	C zustrdr,zustrli,zustrph,
3218	C zvstrdr,zvstrli,zvstrph,
3219	C paprs,u,v,
3220	C aam, torsfc)
3221	ENDIF
3222	cIM cf. FLott END
3223	cIM
3224	IF (ip_ebil_phy.ge.2) THEN
3225	ztit='after orography'
3226	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3227	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3228	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3229	END IF
3230	c
3231	c
3232	!====================================================================
3233	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3234	!====================================================================
3235	! Abderrahmane 24.08.09
3236
3237	IF (ok_cosp) THEN
3238	! adeclarer
3239	#ifdef CPP_COSP
3240	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3241
3242	print*,'freq_cosp',freq_cosp
3243	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3244	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3245	! s ref_liq,ref_ice
3246	call phys_cosp(itap,dtime,freq_cosp,
3247	$ ecrit_mth,ecrit_day,ecrit_hf,overlap,
3248	$ klon,klev,rlon,rlat,presnivs,
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