Context Navigation

physiq.F @ 1425

Last change on this file since 1425 was 1425, checked in by lguez, 14 years ago

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

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

File size: 122.9 KB

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

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

Download in other formats:

Original Format