Context Navigation

physiq_mod.F90 @ 4534

Last change on this file since 4534 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

mise a jour des sources pour ISCCP, CALIPSO et PARASOL
prise en compte des changements de phases pour les nuages (Calipso)
rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

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: 144.2 KB

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

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

Download in other formats:

Original Format