Context Navigation

physiq_mod.F90 @ 2545

Last change on this file since 2545 was 2536, checked in by oboucher, 8 years ago
Introducing stratomask diagnostic on where the stratosphere is flag_aerosol_strat = 2 for CMIP6 strat aerosol forcing
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: 148.7 KB

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

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

Download in other formats:

Original Format