Context Navigation

physiq.F90 @ 2056

Last change on this file since 2056 was 2056, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1997:2055 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 123.6 KB

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

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

Download in other formats:

Original Format