Context Navigation

physiq.F @ 1534

Last change on this file since 1534 was 1534, checked in by musat, 13 years ago

Ajouts CFMIP2/CMIP5

6eme fichier de sortie "stations" histstn.nc qui necessite 2 fichiers: PARAM/npCFMIP_param.data contenant le nombre de points (120 pour simulations AMIP, 73 pour aqua) PARAM/pointlocations.txt contenat le numero, les coordonnees (lon,lat) et le nom de chaque station
flag LOGICAL dans tous les appels histwrite_phy pour pouvoir sortir le fichier histstn.nc

NB: 1) les flags de type phys_ que l'on met dans le physiq.def_L* pour ajouter plus de sorties

necessitent dorenavant 6 valeurs, la 6eme correspondant au fichier histstn.nc

2) par defaut le fichier histstn.nc ne sort pas; pour le sortir ajouter les lignes suivantes

dans physiq.def_L*

### Type de fichier : global (n) ou stations (y)
phys_out_filestations = n n n n n y

introduction de 2 jeux de flags pour les taux des GES; taux actuels avec suffixes _act, taux futurs avec "_per" avec 2 appels au rayonnement si taux "_per" different des taux "_act" (utiles pour diags. CFMIP 4CO2)
flags "betaCRF" pour calculs CRF pour experiences sensibilite proprietes optiques eau liquide nuageuse avec initialisations par defaut; sinon besoin de fichier beta_crf.data

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

File size: 126.9 KB

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

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

Download in other formats:

Original Format