Context Navigation

physiq_mod.F @ 2937

Last change on this file since 2937 was 2836, checked in by abierjon, 3 years ago

VENUS GCM:

INCLUDING THE IONOSPHERE CODE IN VENUS GCM

ATTENTION: INCREASED MEMORY DEMAND

NEEDS AT LEAST 135 GB of allocated memory

==============================================================
===> LIST OF APPENDED FILES AND INTERNAL ADDITIONS <==========
==============================================================

NEW MODEL SPECIES

deftank/traceur-chemistry-IONOSPHERE.def

Neutrals: 4 Species: N, NO, NO2, N(2D)

Ions: 15 species:

CO2+, CO+, O+, O2+, H+,

N2+, H2O+, OH+, C+, HCO+,

H3O+, HCO2+, N+, NO+, elec

FROM 36 to 55 chemical species

NEW KEYWORD OF PHYSIQ.DEF

deftank/physiq-96x96x90-chemistry-IONOSPHERE.def

ok_ionchem: keyword supposed to activate ion chemistry.

(be careful that n, no, no2, n2d and the ion species are in the deftank/traceur-chemistry-IONOSPHERE.def)

ok_jonline: keyword supposed to activate the online photochemistry

(be careful that n, no, no2 and n2d are in the traceur-chemistry-IONOSPHERE.def)

==============================================================
===> LIST OF MODIFIED PROGRAMS <=========================
==============================================================

nonoro_gwd_ran_mod.F90

Change EPFLUXMAX value from 5.E-3 to 1.E-3

photochemistry_venus.F90 (krates; photolysis_online; indices;)

Import of keywords: ok_ionchem, tuneupperatm & ok_jonline
addition of ion species in order
Forcing electroneutrality
Update of the reactions a001 and a002 with taking into account the other species
Change of formula for a002 with the formulation of Baulch et al., 1976 (confirmed by Smith and Robertson, 2008)

photolysis_mod.F90

Modification of the rdsolarflux subroutine to include interpolation with Atlas1 and Atlas3 in connection with E10.7

concentration2.F90

Added chemical species n2d, no, no2 and n in the conductivity calculation.

The ions have been excluded because their sum is 10^{5 times less dense than the neutrals and

their thermal conductivity is unknown}

iono_h.F90

Addition of the phdisrate routine
replace the electronic temperature of Mars by that of

origin = 1: Theis et al. 1980 (Venus) with bilinear interpolation altitude/cos(SZA)
origin = 2: Theis et al. 1984 (Venus) with the formula of the electronic temperature at high solar activity

addition of an ion temperature model based on VIRA

cleshphys.h

added ok_ionchem & ok_jonline in COMMON/clesphys_l/

conf_phys.f90

add ok_ionchem & ok_jonline as parameters to read from physiq.def file set to .false. by default

chemparam_mod.F90

Add species in M_tr and corresponding i_X. Set all i_X to zero before reading traceur-chemistry-IONOSPHERE.def
Added Type_tr table to differentiate species: 1 == neutral, 2 == ION, 3 == liquid, 10 == others

euvheat.F90; hrtherm.F; jthermcalc_e107.F; param_read_e107.F

Normalization with Mars

A.M

File size: 71.2 KB

Line
1	!
2	! $Id: $
3	!
4	MODULE physiq_mod
5
6	IMPLICIT NONE
7
8	CONTAINS
9
10	SUBROUTINE physiq (nlon,nlev,nqmax,
11	. debut,lafin,rjourvrai,gmtime,pdtphys,
12	. paprs,pplay,ppk,pphi,pphis,presnivs,
13	. u,v,t,qx,
14	. flxmw,
15	. d_u, d_v, d_t, d_qx, d_ps)
16
17	c======================================================================
18	c
19	c Modifications pour la physique de Venus
20	c S. Lebonnois (LMD/CNRS) Septembre 2005
21	c
22	c ---------------------------------------------------------------------
23	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
24	c
25	c Objet: Moniteur general de la physique du modele
26	cAA Modifications quant aux traceurs :
27	cAA - uniformisation des parametrisations ds phytrac
28	cAA - stockage des moyennes des champs necessaires
29	cAA en mode traceur off-line
30	c modif ( P. Le Van , 12/10/98 )
31	c
32	c Arguments:
33	c
34	c nlon----input-I-nombre de points horizontaux
35	c nlev----input-I-nombre de couches verticales
36	c nqmax---input-I-nombre de traceurs
37	c debut---input-L-variable logique indiquant le premier passage
38	c lafin---input-L-variable logique indiquant le dernier passage
39	c rjour---input-R-numero du jour de l'experience
40	c gmtime--input-R-fraction de la journee (0 a 1)
41	c pdtphys-input-R-pas d'integration pour la physique (seconde)
42	c paprs---input-R-pression pour chaque inter-couche (en Pa)
43	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
44	c ppk ---input-R-fonction d'Exner au milieu de couche
45	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
46	c pphis---input-R-geopotentiel du sol
47	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
48	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
49	c v-------input-R-vitesse Y (de S a N) en m/s
50	c t-------input-R-temperature (K)
51	c qx------input-R-mass mixing ratio traceurs (kg/kg)
52	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
53	c flxmw---input-R-flux de masse vertical en kg/s
54	c
55	c d_u-----output-R-tendance physique de "u" (m/s/s)
56	c d_v-----output-R-tendance physique de "v" (m/s/s)
57	c d_t-----output-R-tendance physique de "t" (K/s)
58	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
59	c d_ps----output-R-tendance physique de la pression au sol
60	c======================================================================
61	USE ioipsl
62	! USE histcom ! not needed; histcom is included in ioipsl
63	use dimphy
64	USE geometry_mod,only: longitude, latitude, ! in radians
65	& longitude_deg,latitude_deg, ! in degrees
66	& cell_area,dx,dy
67	USE phys_state_var_mod ! Variables sauvegardees de la physique
68	USE cpdet_phy_mod, only: cpdet, t2tpot
69	USE chemparam_mod
70	USE conc
71	USE param_v4_h
72	USE compo_hedin83_mod2
73	use radlwsw_newtoncool_mod, only: radlwsw_newtoncool
74	! use ieee_arithmetic
75	use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy
76	use mod_grid_phy_lmdz, only: nbp_lon
77	use infotrac_phy, only: iflag_trac, tname, ttext
78	use vertical_layers_mod, only: pseudoalt
79	use turb_mod, only : sens, turb_resolved
80	use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
81	use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation
82	use iono_h, only: temp_elect, temp_ion
83	#ifdef CPP_XIOS
84	use xios_output_mod, only: initialize_xios_output,
85	& update_xios_timestep,
86	& send_xios_field
87	use wxios, only: wxios_context_init, xios_context_finalize
88	#endif
89	#ifdef MESOSCALE
90	use comm_wrf
91	#else
92	use iophy
93	use write_field_phy
94	use mod_phys_lmdz_omp_data, ONLY: is_omp_master
95	USE mod_phys_lmdz_para, only : is_parallel,jj_nb,
96	& is_north_pole_phy,
97	& is_south_pole_phy,
98	& is_master
99	#endif
100	IMPLICIT none
101	c======================================================================
102	c CLEFS CPP POUR LES IO
103	c =====================
104	#ifndef MESOSCALE
105	c#define histhf
106	#define histday
107	#define histmth
108	#define histins
109	#endif
110	c======================================================================
111	#include "dimsoil.h"
112	#include "clesphys.h"
113	#include "iniprint.h"
114	#include "timerad.h"
115	#include "tabcontrol.h"
116	#include "nirdata.h"
117	#include "hedin.h"
118	c======================================================================
119	LOGICAL ok_journe ! sortir le fichier journalier
120	save ok_journe
121	c PARAMETER (ok_journe=.true.)
122	c
123	LOGICAL ok_mensuel ! sortir le fichier mensuel
124	save ok_mensuel
125	c PARAMETER (ok_mensuel=.true.)
126	c
127	LOGICAL ok_instan ! sortir le fichier instantane
128	save ok_instan
129	c PARAMETER (ok_instan=.true.)
130	c
131	c======================================================================
132	c
133	c Variables argument:
134	c
135	INTEGER nlon
136	INTEGER nlev
137	INTEGER nqmax
138	REAL rjourvrai
139	REAL gmtime
140	REAL pdtphys
141	LOGICAL debut, lafin
142	REAL paprs(klon,klev+1)
143	REAL pplay(klon,klev)
144	REAL pphi(klon,klev)
145	REAL pphis(klon)
146	REAL presnivs(klev)
147
148	! ADAPTATION GCM POUR CP(T)
149	REAL ppk(klon,klev)
150
151	REAL u(klon,klev)
152	REAL v(klon,klev)
153	REAL t(klon,klev)
154	REAL qx(klon,klev,nqmax)
155
156	REAL d_u_dyn(klon,klev)
157	REAL d_t_dyn(klon,klev)
158
159	REAL flxmw(klon,klev)
160
161	REAL d_u(klon,klev)
162	REAL d_v(klon,klev)
163	REAL d_t(klon,klev)
164	REAL d_qx(klon,klev,nqmax)
165	REAL d_ps(klon)
166
167	logical ok_hf
168	real ecrit_hf
169	integer nid_hf
170	save ok_hf, ecrit_hf, nid_hf
171
172	#ifdef histhf
173	data ok_hf,ecrit_hf/.true.,0.25/
174	#else
175	data ok_hf/.false./
176	#endif
177
178	c Variables propres a la physique
179
180	integer,save :: itap ! physics counter
181	REAL delp(klon,klev) ! epaisseur d'une couche
182	REAL omega(klon,klev) ! vitesse verticale en Pa/s (+ downward)
183	REAL vertwind(klon,klev) ! vitesse verticale en m/s (+ upward)
184
185	INTEGER igwd,idx(klon),itest(klon)
186
187	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
188
189	REAL zulow(klon),zvlow(klon)
190	REAL zustrdr(klon), zvstrdr(klon)
191	REAL zustrli(klon), zvstrli(klon)
192	REAL zustrhi(klon), zvstrhi(klon)
193	REAL zublstrdr(klon), zvblstrdr(klon)
194	REAL znlow(klon), zeff(klon)
195	REAL zbl(klon), knu2(klon),kbreak(nlon)
196	REAL tau0(klon), ztau(klon,klev)
197
198	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
199	real zdtlev(klon,klev),zdzlev(klon,klev)
200	real ztlev(klon,klev)
201	real zn2(klon,klev) ! BV^2 at plev
202
203	c Pour les bilans de moment angulaire,
204	integer bilansmc
205	c Pour le transport de ballons
206	integer ballons
207	c j'ai aussi besoin
208	c du stress de couche limite a la surface:
209
210	REAL zustrcl(klon),zvstrcl(klon)
211
212	c et du stress total c de la physique:
213
214	REAL zustrph(klon),zvstrph(klon)
215
216	c Variables locales:
217	c
218	REAL cdragh(klon) ! drag coefficient pour T and Q
219	REAL cdragm(klon) ! drag coefficient pour vent
220	c
221	cAA Pour TRACEURS
222	cAA
223	REAL,save,allocatable :: source(:,:)
224	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
225	REAL yu1(klon) ! vents dans la premiere couche U
226	REAL yv1(klon) ! vents dans la premiere couche V
227
228	REAL dsens(klon) ! derivee chaleur sensible
229	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
230	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
231	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
232	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
233	c
234	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
235
236	c======================================================================
237	c
238	c Declaration des procedures appelees
239	c
240	EXTERNAL ajsec ! ajustement sec
241	EXTERNAL clmain ! couche limite
242	EXTERNAL clmain_ideal ! couche limite simple
243	EXTERNAL hgardfou ! verifier les temperatures
244	c EXTERNAL orbite ! calculer l'orbite
245	EXTERNAL phyetat0 ! lire l'etat initial de la physique
246	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
247	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
248	! EXTERNAL suphec ! initialiser certaines constantes
249	EXTERNAL transp ! transport total de l'eau et de l'energie
250	EXTERNAL printflag
251	EXTERNAL zenang
252	EXTERNAL diagetpq
253	EXTERNAL conf_phys
254	EXTERNAL diagphy
255	EXTERNAL mucorr
256	EXTERNAL nirco2abs
257	EXTERNAL nir_leedat
258	EXTERNAL nltecool
259	EXTERNAL nlte_tcool
260	EXTERNAL nlte_setup
261	EXTERNAL blendrad
262	EXTERNAL nlthermeq
263	EXTERNAL euvheat
264	EXTERNAL param_read_e107
265	EXTERNAL conduction
266	EXTERNAL molvis
267	EXTERNAL moldiff_red
268
269	c
270	c Variables locales
271	c
272	CXXX PB
273	REAL fluxt(klon,klev) ! flux turbulent de chaleur
274	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
275	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
276	c
277	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
278	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
279	REAL flux_ec(klon,klev) ! flux de chaleur Ec
280	c
281	REAL tmpout(klon,klev) ! [K/s]
282	c
283	REAL dist, rmu0(klon), fract(klon)
284	REAL zdtime, zlongi
285	c
286	INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev, isoil
287	c
288	REAL zphi(klon,klev)
289	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
290	REAL tlaymean ! valeur temporaire pour calculer zzlay
291	real tsurf(klon)
292
293	c va avec nlte_tcool
294	INTEGER ierr_nlte
295	REAL varerr
296
297	! photochemistry
298
299	integer :: chempas
300	real :: nbapp_chem, zctime
301
302	! sedimentation
303
304	REAL :: m0_mode1(klev,2),m0_mode2(klev,2)
305	REAL :: m3_mode1(klev,3),m3_mode2 (klev,3)
306	REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2)
307	REAL :: aer_flux(klev)
308	c
309	c Variables du changement
310	c
311	c ajs: ajustement sec
312	c vdf: couche limite (Vertical DiFfusion)
313	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
314	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
315	c
316	REAL d_ts(klon)
317	c
318	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
319	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
320	c
321	CMOD LOTT: Tendances Orography Sous-maille
322	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
323	REAL d_t_oro(klon,klev)
324	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
325	REAL d_t_lif(klon,klev)
326	C Tendances Ondes de G non oro (runs strato).
327	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
328	REAL d_t_hin(klon,klev)
329
330	c Tendencies due to radiative scheme [K/s]
331	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
332	c are not computed at each physical timestep
333	c therefore, they are defined and saved in phys_state_var_mod
334
335	c Tendencies due to molecular viscosity and conduction
336	real d_t_conduc(klon,klev) ! [K/s]
337	real d_u_molvis(klon,klev) ! [m/s] /s
338	real d_v_molvis(klon,klev) ! [m/s] /s
339
340	c Tendencies due to molecular diffusion
341	real d_q_moldif(klon,klev,nqmax)
342
343	c
344	c Variables liees a l'ecriture de la bande histoire physique
345	c
346	INTEGER ecrit_mth
347	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
348	c
349	INTEGER ecrit_day
350	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
351	c
352	INTEGER ecrit_ins
353	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
354	c
355	integer itau_w ! pas de temps ecriture = itap + itau_phy
356
357	c Variables locales pour effectuer les appels en serie
358	c
359	REAL t_seri(klon,klev)
360	REAL u_seri(klon,klev), v_seri(klon,klev)
361	c
362	REAL :: tr_seri(klon,klev,nqmax)
363	REAL :: tr_hedin(klon,klev,nqmax)
364	REAL :: d_tr(klon,klev,nqmax)
365	c pour sorties
366	REAL :: col_dens_tr(klon,nqmax)
367	REAL,allocatable,save :: prod_tr(:,:,:)
368	REAL,allocatable,save :: loss_tr(:,:,:)
369
370	c pour ioipsl
371	INTEGER nid_day, nid_mth, nid_ins
372	SAVE nid_day, nid_mth, nid_ins
373	INTEGER nhori, nvert, idayref
374	REAL zsto, zout, zsto1, zsto2, zero
375	parameter (zero=0.0e0)
376	real zjulian
377	save zjulian
378
379	CHARACTER*2 str2
380	character*20 modname
381	character*80 abort_message
382	logical ok_sync
383
384	character*30 nom_fichier
385	character*10 varname
386	character*40 vartitle
387	character*20 varunits
388	C Variables liees au bilan d'energie et d'enthalpi
389	REAL ztsol(klon)
390	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
391	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
392	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
393	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
394	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
395	REAL d_h_vcol_phy
396	REAL fs_bound, fq_bound
397	SAVE d_h_vcol_phy
398	REAL zero_v(klon),zero_v2(klon,klev)
399	CHARACTER*15 ztit
400	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
401	SAVE ip_ebil
402	DATA ip_ebil/2/
403	INTEGER if_ebil ! level for energy conserv. dignostics
404	SAVE if_ebil
405	c+jld ec_conser
406	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
407	c-jld ec_conser
408
409	c ALBEDO VARIATIONS (VCD)
410	REAL factAlb
411	c TEST VENUS...
412	REAL mang(klon,klev) ! moment cinetique
413	REAL mangtot ! moment cinetique total
414
415	c cell_area for outputs in hist*
416	REAL cell_area_out(klon)
417	#ifdef MESOSCALE
418	REAL :: dt_dyn(klev)
419	#endif
420	c Declaration des constantes et des fonctions thermodynamiques
421	c
422	#include "YOMCST.h"
423
424	c======================================================================
425	c======================================================================
426	c INITIALISATIONS
427	c======================================================================
428
429	modname = 'physiq'
430	ok_sync=.TRUE.
431
432	bilansmc = 0
433	ballons = 0
434	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
435	#ifndef MESOSCALE
436	if (is_parallel) then
437	bilansmc = 0
438	ballons = 0
439	endif
440	#endif
441	IF (if_ebil.ge.1) THEN
442	DO i=1,klon
443	zero_v(i)=0.
444	END DO
445	DO i=1,klon
446	DO j=1,klev
447	zero_v2(i,j)=0.
448	END DO
449	END DO
450	END IF
451
452	c======================================================================
453	c PREMIER APPEL SEULEMENT
454	c========================
455	IF (debut) THEN
456	allocate(source(klon,nqmax))
457	allocate(prod_tr(klon,klev,nqmax))
458	allocate(loss_tr(klon,klev,nqmax))
459	allocate(no_emission(klon,klev))
460	allocate(o2_emission(klon,klev))
461
462	#ifdef CPP_XIOS
463	! Initialize XIOS context
464	write(,) "physiq: call wxios_context_init"
465	CALL wxios_context_init
466	#endif
467
468	! The call to suphec is now done in iniphysiq_mod (interface)
469	! CALL suphec ! initialiser constantes et parametres phys.
470
471	IF (if_ebil.ge.1) d_h_vcol_phy=0.
472	c
473	c appel a la lecture du physiq.def
474	c
475	call conf_phys(ok_journe, ok_mensuel,
476	. ok_instan,
477	. if_ebil)
478
479	call phys_state_var_init(nqmax)
480	c
481	c Initialising Hedin model for upper atm
482	c (to be revised when coupled to chemistry) :
483	call conc_init
484
485	! initialise physics counter
486
487	itap = 0
488
489	#ifdef MESOSCALE
490	print*,'check pdtphys',pdtphys
491	PRINT*,'check phisfi ',pphis(1),pphis(klon)
492	PRINT*,'check geop',pphi(1,1),pphi(klon,klev)
493	PRINT*,'check radsol',radsol(1),radsol(klon)
494	print*,'check ppk',ppk(1,1),ppk(klon,klev)
495	print*,'check ftsoil',ftsoil(1,1),ftsoil(klon,nsoilmx)
496	print*,'check ftsol',ftsol(1),ftsol(klon)
497	print*, "check temp", t(1,1),t(klon,klev)
498	print*, "check pres",paprs(1,1),paprs(klon,klev),pplay(1,1),
499	. pplay(klon,klev)
500	print*, "check u", u(1,1),u(klon,klev)
501	print*, "check v", v(1,1),v(klon,klev)
502	print*,'check falbe',falbe(1),falbe(klon)
503	!nqtot=nqmax
504	!ALLOCATE(tname(nqtot))
505	!tname=noms
506	zmea=0.
507	zstd=0.
508	zsig=0.
509	zgam=0.
510	zthe=0.
511	dtime=pdtphys
512	#else
513	c
514	c Lecture startphy.nc :
515	c
516	CALL phyetat0 ("startphy.nc")
517	IF (.not.startphy_file) THEN
518	! Additionnal academic initializations
519	ftsol(:)=t(:,1) ! surface temperature as in first atm. layer
520	DO isoil=1, nsoilmx
521	! subsurface temperatures equal to surface temperature
522	ftsoil(:,isoil)=ftsol(:)
523	ENDDO
524	ENDIF
525	#endif
526
527	c dtime est defini dans tabcontrol.h et lu dans startphy
528	c pdtphys est calcule a partir des nouvelles conditions:
529	c Reinitialisation du pas de temps physique quand changement
530	IF (ABS(dtime-pdtphys).GT.0.001) THEN
531	WRITE(lunout,*) 'Pas physique a change',dtime,
532	. pdtphys
533	c abort_message='Pas physique n est pas correct '
534	c call abort_physic(modname,abort_message,1)
535	c----------------
536	c pour initialiser convenablement le time_counter, il faut tenir compte
537	c du changement de dtime en changeant itau_phy (point de depart)
538	itau_phy = NINT(itau_phy*dtime/pdtphys)
539	c----------------
540	dtime=pdtphys
541	ENDIF
542
543	radpas = NINT(RDAY/pdtphys/nbapp_rad)
544
545	CALL printflag( ok_journe,ok_instan )
546
547	#ifdef CPP_XIOS
548	write(,) "physiq: call initialize_xios_output"
549	call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY,
550	& presnivs,pseudoalt)
551	#endif
552
553	c
554	c---------
555	c FLOTT
556	IF (ok_orodr) THEN
557	DO i=1,klon
558	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
559	ENDDO
560	CALL SUGWD(klon,klev,paprs,pplay)
561	DO i=1,klon
562	zuthe(i)=0.
563	zvthe(i)=0.
564	if(zstd(i).gt.10.)then
565	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
566	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
567	endif
568	ENDDO
569	ENDIF
570
571	if (bilansmc.eq.1) then
572	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
573	C DU BILAN DE MOMENT ANGULAIRE.
574	open(27,file='aaam_bud.out',form='formatted')
575	open(28,file='fields_2d.out',form='formatted')
576	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
577	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
578	endif !bilansmc
579
580	c--------------SLEBONNOIS
581	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
582	C DES BALLONS
583	if (ballons.eq.1) then
584	open(30,file='ballons-lat.out',form='formatted')
585	open(31,file='ballons-lon.out',form='formatted')
586	open(32,file='ballons-u.out',form='formatted')
587	open(33,file='ballons-v.out',form='formatted')
588	open(34,file='ballons-alt.out',form='formatted')
589	write(,)'Ouverture des ballons*.out'
590	endif !ballons
591	c-------------
592
593	c---------
594	C TRACEURS
595	C source dans couche limite
596	source(:,:) = 0.0 ! pas de source, pour l'instant
597	c---------
598
599	c---------
600	c INITIALIZE THERMOSPHERIC PARAMETERS
601	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
602
603	if (callthermos) then
604	call fill_data_thermos
605	call allocate_param_thermos(klev)
606	call param_read_e107
607	endif
608
609	c Initialisation (recomputed in concentration2)
610	do ig=1,klon
611	do j=1,klev
612	rnew(ig,j)=RD
613	cpnew(ig,j)=cpdet(t(ig,j))
614	mmean(ig,j)=RMD
615	akknew(ig,j)=1.e-4
616	rho(ig,j)=pplay(ig,j)/(rnew(ig,j)*t(ig,j))
617	enddo
618
619	enddo
620
621	IF(callthermos.or.callnlte.or.callnirco2) THEN
622	call compo_hedin83_init2
623	ENDIF
624	if (callnlte.and.nltemodel.eq.2) call nlte_setup
625	if (callnirco2.and.nircorr.eq.1) call nir_leedat
626	c---------
627
628	c
629	c Verifications:
630	c
631	IF (nlon .NE. klon) THEN
632	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
633	. klon
634	abort_message='nlon et klon ne sont pas coherents'
635	call abort_physic(modname,abort_message,1)
636	ENDIF
637	IF (nlev .NE. klev) THEN
638	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
639	. klev
640	abort_message='nlev et klev ne sont pas coherents'
641	call abort_physic(modname,abort_message,1)
642	ENDIF
643	c
644	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
645	$ THEN
646	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
647	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
648	abort_message='Nbre d appels au rayonnement insuffisant'
649	call abort_physic(modname,abort_message,1)
650	ENDIF
651	c
652	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
653	. iflag_ajs
654	c
655	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
656	IF (ok_mensuel) THEN
657	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
658	. ecrit_mth
659	ENDIF
660
661	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
662	IF (ok_journe) THEN
663	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
664	. ecrit_day
665	ENDIF
666
667	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
668	IF (ok_instan) THEN
669	WRITE(lunout,*)'La frequence de sortie instant. est de ',
670	. ecrit_ins
671	ENDIF
672
673	c Initialisation des sorties
674	c===========================
675
676	#ifdef CPP_IOIPSL
677
678	#ifdef histhf
679	#include "ini_histhf.h"
680	#endif
681
682	#ifdef histday
683	#include "ini_histday.h"
684	#endif
685
686	#ifdef histmth
687	#include "ini_histmth.h"
688	#endif
689
690	#ifdef histins
691	#include "ini_histins.h"
692	#endif
693
694	#endif
695
696	c
697	c Initialiser les valeurs de u pour calculs tendances
698	c (pour T, c'est fait dans phyetat0)
699	c
700	DO k = 1, klev
701	DO i = 1, klon
702	u_ancien(i,k) = u(i,k)
703	ENDDO
704	ENDDO
705
706	! initialisation of microphysical and chemical parameters
707
708	if (ok_chem .and. .not. ok_cloud) then
709	print*,"chemistry requires clouds"
710	print*,"ok_cloud must be .true."
711	stop
712	end if
713
714	if (.not. ok_chem .and. ok_cloud .and. (cl_scheme == 1)) then
715	print*,"cl_scheme = 1 requires chemistry"
716	print*,"ok_chem must be .true."
717	stop
718	end if
719
720	! number of microphysical tracers
721
722	nmicro = 0
723	if (ok_cloud .and. (cl_scheme == 1)) nmicro = 2
724	if (ok_cloud .and. (cl_scheme == 2)) nmicro = 12
725
726	! initialise chemical parameters. includes the indexation of microphys tracers
727
728	if (ok_chem .or. cl_scheme == 2) then
729	call chemparam_ini()
730	end if
731
732	! initialise cloud parameters (for cl_scheme = 1)
733
734	if (ok_cloud .and. (cl_scheme == 1)) then
735	call cloud_ini(nlon,nlev)
736	end if
737
738	! initialise mmean
739
740	if (callthermos) then
741	call concentrations2(pplay,t,qx,nqmax)
742	end if
743
744	c========================
745	ENDIF ! debut
746	c========================
747
748	c======================================================================
749	! ------------------------------------------------------
750	! Initializations done at every physical timestep:
751	! ------------------------------------------------------
752
753	c Mettre a zero des variables de sortie (pour securite)
754	c
755	DO i = 1, klon
756	d_ps(i) = 0.0
757	ENDDO
758	DO k = 1, klev
759	DO i = 1, klon
760	d_t(i,k) = 0.0
761	d_u(i,k) = 0.0
762	d_v(i,k) = 0.0
763	ENDDO
764	ENDDO
765	DO iq = 1, nqmax
766	DO k = 1, klev
767	DO i = 1, klon
768	d_qx(i,k,iq) = 0.0
769	ENDDO
770	ENDDO
771	ENDDO
772	c
773	c Ne pas affecter les valeurs entrees de u, v, h, et q
774	c
775	DO k = 1, klev
776	DO i = 1, klon
777	t_seri(i,k) = t(i,k)
778	u_seri(i,k) = u(i,k)
779	v_seri(i,k) = v(i,k)
780	ENDDO
781	ENDDO
782	DO iq = 1, nqmax
783	DO k = 1, klev
784	DO i = 1, klon
785	tr_seri(i,k,iq) = qx(i,k,iq)
786	ENDDO
787	ENDDO
788	ENDDO
789	C
790	DO i = 1, klon
791	ztsol(i) = ftsol(i)
792	ENDDO
793	C
794	IF (if_ebil.ge.1) THEN
795	ztit='after dynamic'
796	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
797	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
798	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
799	C Comme les tendances de la physique sont ajoute dans la dynamique,
800	C on devrait avoir que la variation d'entalpie par la dynamique
801	C est egale a la variation de la physique au pas de temps precedent.
802	C Donc la somme de ces 2 variations devrait etre nulle.
803	call diagphy(cell_area,ztit,ip_ebil
804	e , zero_v, zero_v, zero_v, zero_v, zero_v
805	e , zero_v, zero_v, zero_v, ztsol
806	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
807	s , fs_bound, fq_bound )
808	END IF
809
810	c====================================================================
811	c XIOS outputs
812
813	#ifdef CPP_XIOS
814	! update XIOS time/calendar
815	call update_xios_timestep
816	#endif
817
818	c====================================================================
819	c Diagnostiquer la tendance dynamique
820	c
821	IF (ancien_ok) THEN
822	DO k = 1, klev
823	DO i = 1, klon
824	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
825	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
826	ENDDO
827	ENDDO
828
829	! ADAPTATION GCM POUR CP(T)
830	do i=1,klon
831	flux_dyn(i,1) = 0.0
832	do j=2,klev
833	flux_dyn(i,j) = flux_dyn(i,j-1)
834	. +cpnew(i,j-1)/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
835	enddo
836	enddo
837
838	ELSE
839	DO k = 1, klev
840	DO i = 1, klon
841	d_u_dyn(i,k) = 0.0
842	d_t_dyn(i,k) = 0.0
843	ENDDO
844	ENDDO
845	ancien_ok = .TRUE.
846	ENDIF
847	c====================================================================
848
849	! Compute vertical velocity (Pa/s) from vertical mass flux
850	! Need to linearly interpolate mass flux to mid-layers
851	do k=1,klev-1
852	omega(1:klon,k) = 0.5RG(flxmw(1:klon,k)+flxmw(1:klon,k+1))
853	. / cell_area(1:klon)
854	enddo
855	omega(1:klon,klev) = 0.5RGflxmw(1:klon,klev) / cell_area(1:klon)
856
857	c======
858	c GEOP CORRECTION
859	c
860	c Ajouter le geopotentiel du sol:
861	c
862	DO k = 1, klev
863	DO i = 1, klon
864	zphi(i,k) = pphi(i,k) + pphis(i)
865	ENDDO
866	ENDDO
867
868	c............................
869	c CETTE CORRECTION VA DE PAIR AVEC DES MODIFS DE LEAPFROG(_p)
870	c ELLE MARCHE A 50 NIVEAUX (si mmean constante...)
871	c MAIS PAS A 78 NIVEAUX (quand mmean varie...)
872	c A ANALYSER PLUS EN DETAIL AVANT D'UTILISER
873	c............................
874	c zphi is recomputed (pphi is not ok if mean molecular mass varies)
875	c with dphi = RT/mmean d(ln p) [evaluated at interface]
876
877	c DO i = 1, klon
878	c zphi(i,1) = pphis(i) + Rt_seri(i,1)/mmean(i,1)1000.
879	c * *( log(paprs(i,1)) - log(pplay(i,1)) )
880	c DO k = 2, klev
881	c zphi(i,k) = zphi(i,k-1)
882	c * + R500.(t_seri(i,k)/mmean(i,k)+t_seri(i,k-1)/mmean(i,k-1))
883	c * * (log(pplay(i,k-1)) - log(pplay(i,k)))
884	c ENDDO
885	c ENDDO
886	c............................
887	c=====
888
889	c calcul de l'altitude aux niveaux intercouches
890	c ponderation des altitudes au niveau des couches en dp/p
891
892	DO i=1,klon
893	zzlay(i,1)=zphi(i,1)/RG ! [m]
894	zzlev(i,1)=pphis(i)/RG ! [m]
895	ENDDO
896	DO k=2,klev
897	DO i=1,klon
898	tlaymean=t_seri(i,k)
899	IF (t_seri(i,k).ne.t_seri(i,k-1))
900	& tlaymean=(t_seri(i,k)-t_seri(i,k-1))
901	& /log(t_seri(i,k)/t_seri(i,k-1))
902
903	zzlay(i,k)=zzlay(i,k-1)
904	& -(log(pplay(i,k)/pplay(i,k-1))rnew(i,k-1)tlaymean
905	& /(RG(RA/(RA+zzlay(i,k-1)))*2))
906	ENDDO
907	ENDDO
908	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
909	! Old : from geopotential. Problem when mu varies in the upper atm...
910	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
911	! DO k=1,klev
912	! DO i=1,klon
913	! zzlay(i,k)=zphi(i,k)/RG ! [m]
914	! ENDDO
915	! ENDDO
916	! DO i=1,klon
917	! zzlev(i,1)=pphis(i)/RG ! [m]
918	! ENDDO
919	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
920	DO k=2,klev
921	DO i=1,klon
922	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
923	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
924	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
925	ENDDO
926	ENDDO
927	DO i=1,klon
928	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
929	ENDDO
930
931	c====================================================================
932	c
933	c Verifier les temperatures
934	c
935	CALL hgardfou(t_seri,ftsol,'debutphy')
936
937	c====================================================================
938	c Orbite et eclairement
939	c=======================
940
941	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
942	c donc pas de variations de Ls, ni de dist.
943	c La seule chose qui compte, c'est la rotation de la planete devant
944	c le Soleil...
945
946	zlongi = 0.0
947	dist = 0.72333 ! en UA
948
949	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
950	c a sa valeur, et prendre en compte leur evolution,
951	c il faudra refaire un orbite.F...
952	c CALL orbite(zlongi,dist)
953
954	IF (cycle_diurne) THEN
955	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
956	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
957	& rmu0,fract)
958	ELSE
959	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
960	ENDIF
961
962	! print fraction of venus day
963
964	if (is_master) then
965	print*, 'gmtime = ', gmtime
966	end if
967
968	c======================================================================
969	c FIN INITIALISATIONS
970	c======================================================================
971	c======================================================================
972
973	c=======================================================
974	! CONDUCTION and MOLECULAR VISCOSITY
975	c=======================================================
976
977	d_t_conduc(:,:)=0.
978	d_u_molvis(:,:)=0.
979	d_v_molvis(:,:)=0.
980
981	IF (callthermos) THEN
982
983	tsurf(:)=t_seri(:,1)
984	call conduction(klon, klev,pdtphys,
985	$ pplay,paprs,t_seri,
986	$ tsurf,zzlev,zzlay,d_t_conduc)
987
988	call molvis(klon, klev,pdtphys,
989	$ pplay,paprs,t_seri,
990	$ u,tsurf,zzlev,zzlay,d_u_molvis)
991
992	call molvis(klon, klev, pdtphys,
993	$ pplay,paprs,t_seri,
994	$ v,tsurf,zzlev,zzlay,d_v_molvis)
995
996	DO k=1,klev
997	DO ig=1,klon
998	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
999	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! [m/s]
1000	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! [m/s]
1001	ENDDO
1002	ENDDO
1003	ENDIF
1004
1005	c====================================================================
1006	c Compute mean mass, cp and R :
1007	c------------------------------------
1008
1009	if(callthermos) then
1010	call concentrations2(pplay,t_seri,tr_seri, nqmax)
1011	endif
1012
1013
1014	c=======================================================
1015	! CHEMISTRY AND MICROPHYSICS
1016	c=======================================================
1017
1018	if (iflag_trac.eq.1) then
1019	!====================================================================
1020	! Case 1: pseudo-chemistry with relaxation toward fixed profile
1021	!=========
1022	if (tr_scheme.eq.1) then
1023
1024	call phytrac_relax (debut,lafin,nqmax,
1025	I nlon,nlev,dtime,pplay,
1026	O tr_seri)
1027
1028	elseif (tr_scheme.eq.2) then
1029	!====================================================================
1030	! Case 2: surface emission
1031	! For the moment, inspired from Mars version
1032	! However, the variable 'source' could be used in physiq
1033	! so the call to phytrac_emiss could be to initialise it.
1034	!=========
1035
1036	call phytrac_emiss (debut,lafin,nqmax,
1037	I nlon,nlev,dtime,paprs,
1038	I latitude_deg,longitude_deg,
1039	O tr_seri)
1040
1041	else if (tr_scheme.eq.3) then
1042	!====================================================================
1043	! Case 3: Full chemistry and/or clouds.
1044	! routines are called every "chempas" physical timestep.
1045	!
1046	! if the physics is called 96000 times per venus day:
1047	!
1048	! nbapp_chem = 24000 => chempas = 4 => zctime = 420 s
1049	! nbapp_chem = 12000 => chempas = 8 => zctime = 840 s
1050	!=========
1051
1052	nbapp_chem = 24000
1053	chempas = nint(rday/pdtphys/nbapp_chem)
1054	zctime = dtime*real(chempas) ! chemical timestep
1055
1056	if (mod(itap,chempas) == 0) then ! <------- start of chemistry supercycling
1057
1058	! photochemistry and microphysics
1059
1060	call phytrac_chimie(debut,
1061	$ gmtime,
1062	$ nqmax,
1063	$ klon,
1064	$ latitude_deg,
1065	$ longitude_deg,
1066	$ zzlay,
1067	$ nlev,
1068	$ zctime,
1069	$ t_seri,
1070	$ pplay,
1071	$ tr_seri,
1072	$ d_tr_chem,
1073	$ iter,
1074	$ prod_tr,
1075	$ loss_tr,
1076	$ no_emission,
1077	$ o2_emission)
1078
1079	if (ok_sedim) then
1080	if (cl_scheme == 1) then
1081
1082	! sedimentation for simplified microphysics
1083
1084	#ifndef MESOSCALE
1085	call new_cloud_sedim(klon,
1086	$ nlev,
1087	$ zctime,
1088	$ pplay,
1089	$ paprs,
1090	$ t_seri,
1091	$ tr_seri,
1092	$ d_tr_chem,
1093	$ d_tr_sed(:,:,1:2),
1094	$ nqmax,
1095	$ Fsedim)
1096
1097	! test to avoid nans
1098
1099	do k = 1, klev
1100	do i = 1, klon
1101	if ((d_tr_sed(i,k,1) /= d_tr_sed(i,k,1)) .or.
1102	$ (d_tr_sed(i,k,2) /= d_tr_sed(i,k,2))) then
1103	print*,'sedim NaN PROBLEM'
1104	print*,'d_tr_sed Nan?',d_tr_sed(i,k,:)
1105	print*,'Temp',t_seri(i,k)
1106	print*,'lat-lon',i,'level',k,'zctime',zctime
1107	print*,'F_sed',Fsedim(i,k)
1108	d_tr_sed(i,k,:) = 0.
1109	end if
1110	end do
1111	end do
1112
1113	! tendency due to condensation and sedimentation
1114
1115	d_tr_sed(:,:,1:2) = d_tr_sed(:,:,1:2)/zctime
1116	Fsedim(:,1:klev) = Fsedim(:,1:klev)/zctime
1117	Fsedim(:,klev+1) = 0.
1118
1119	else if (cl_scheme == 2) then
1120
1121	! sedimentation for detailed microphysics
1122
1123	d_tr_sed(:,:,:) = 0.
1124
1125	do i = 1, klon
1126
1127	! mode 1
1128
1129	m0_mode1(:,1) = tr_seri(i,:,i_m0_mode1drop)
1130	m0_mode1(:,2) = tr_seri(i,:,i_m0_mode1ccn)
1131	m3_mode1(:,1) = tr_seri(i,:,i_m3_mode1sa)
1132	m3_mode1(:,2) = tr_seri(i,:,i_m3_mode1w)
1133	m3_mode1(:,3) = tr_seri(i,:,i_m3_mode1ccn)
1134
1135	call drop_sedimentation(zctime,klev,m0_mode1,m3_mode1,
1136	$ paprs(i,:),zzlev(i,:),
1137	$ zzlay(i,:),t_seri(i,:),1,
1138	$ d_ccn_sed(:,1),d_drop_sed,
1139	$ d_ccn_sed(:,2),d_liq_sed)
1140
1141	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
1142	$ + d_drop_sed
1143	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
1144	$ + d_ccn_sed(:,1)
1145	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
1146	$ + d_ccn_sed(:,2)
1147	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
1148	$ + d_liq_sed(:,1)
1149	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
1150	$ + d_liq_sed(:,2)
1151
1152	! mode 2
1153
1154	m0_mode2(:,1) = tr_seri(i,:,i_m0_mode2drop)
1155	m0_mode2(:,2) = tr_seri(i,:,i_m0_mode2ccn)
1156	m3_mode2(:,1) = tr_seri(i,:,i_m3_mode2sa)
1157	m3_mode2(:,2) = tr_seri(i,:,i_m3_mode2w)
1158	m3_mode2(:,3) = tr_seri(i,:,i_m3_mode2ccn)
1159
1160	call drop_sedimentation(zctime,klev,m0_mode2,m3_mode2,
1161	$ paprs(i,:),zzlev(i,:),
1162	$ zzlay(i,:),t_seri(i,:),2,
1163	$ d_ccn_sed(:,1),d_drop_sed,
1164	$ d_ccn_sed(:,2),d_liq_sed)
1165
1166	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
1167	$ + d_drop_sed
1168	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
1169	$ + d_ccn_sed(:,1)
1170	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
1171	$ + d_ccn_sed(:,2)
1172	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
1173	$ + d_liq_sed(:,1)
1174	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
1175	$ + d_liq_sed(:,2)
1176
1177	! aer
1178
1179	call aer_sedimentation(zctime,klev,
1180	$ tr_seri(i,:,i_m0_aer),
1181	$ tr_seri(i,:,i_m3_aer),
1182	$ paprs(i,:),zzlev(i,:),
1183	$ zzlay(i,:),t_seri(i,:),
1184	$ d_tr_sed(i,:,i_m0_aer),
1185	$ d_tr_sed(i,:,i_m3_aer),
1186	$ aer_flux)
1187
1188	end do
1189
1190	! tendency due to sedimentation
1191
1192	do iq = nqmax-nmicro+1,nqmax
1193	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime
1194	end do
1195	#endif
1196	end if ! cl_scheme
1197
1198	! update gaseous tracers (chemistry)
1199
1200	do iq = 1, nqmax - nmicro
1201	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
1202	$ + d_tr_chem(:,:,iq)*zctime,1.e-30)
1203	end do
1204
1205	! update condensed tracers (condensation + sedimentation)
1206
1207	if (cl_scheme == 1) then
1208	tr_seri(:,:,i_h2so4liq) = max(tr_seri(:,:,i_h2so4liq)
1209	$ + d_tr_sed(:,:,1)*zctime, 1.e-30)
1210	tr_seri(:,:,i_h2oliq) = max(tr_seri(:,:,i_h2oliq)
1211	$ + d_tr_sed(:,:,2)*zctime, 1.e-30)
1212	else if (cl_scheme == 2) then
1213	do iq = nqmax-nmicro+1,nqmax
1214	tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
1215	$ + d_tr_sed(:,:,iq)*zctime,1.e-30)
1216	end do
1217	end if ! cl_scheme
1218
1219	end if ! ok_sedim
1220	end if ! mod(itap,chempas) <------- end of chemistry supercycling
1221
1222	!=========
1223	! End Case 3: Full chemistry and/or clouds.
1224	!====================================================================
1225
1226	end if ! tr_scheme
1227	end if ! iflag_trac
1228
1229	c====================================================================
1230	c Appeler la diffusion verticale (programme de couche limite)
1231	c====================================================================
1232
1233	c-------------------------------
1234	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
1235	c l'equilibre radiatif du sol
1236	if (.not. ok_clmain) then
1237	if (debut) then
1238	print*,"ATTENTION, CLMAIN SHUNTEE..."
1239	endif
1240
1241	DO i = 1, klon
1242	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
1243	fder(i) = 0.0e0
1244	dlw(i) = 0.0e0
1245	ENDDO
1246
1247	c Incrementer la temperature du sol
1248	c
1249	DO i = 1, klon
1250	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
1251	ftsol(i) = ftsol(i) + d_ts(i)
1252	do j=1,nsoilmx
1253	ftsoil(i,j)=ftsol(i)
1254	enddo
1255	ENDDO
1256
1257	c-------------------------------
1258	else
1259	c-------------------------------
1260
1261	fder = dlw
1262
1263	! ADAPTATION GCM POUR CP(T)
1264
1265	if (physideal) then
1266	CALL clmain_ideal(dtime,itap,
1267	e t_seri,u_seri,v_seri,
1268	e rmu0,
1269	e ftsol,
1270	$ ftsoil,
1271	$ paprs,pplay,ppk,radsol,falbe,
1272	e solsw, sollw, sollwdown, fder,
1273	e longitude_deg, latitude_deg, dx, dy,
1274	e debut, lafin,
1275	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1276	s fluxt,fluxu,fluxv,cdragh,cdragm,
1277	s dsens,
1278	s ycoefh,yu1,yv1)
1279	else
1280	CALL clmain(dtime,itap,
1281	e t_seri,u_seri,v_seri,
1282	e rmu0,
1283	e ftsol,
1284	$ ftsoil,
1285	$ paprs,pplay,ppk,radsol,falbe,
1286	e solsw, sollw, sollwdown, fder,
1287	e longitude_deg, latitude_deg, dx, dy,
1288	& q2,
1289	e debut, lafin,
1290	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1291	s fluxt,fluxu,fluxv,cdragh,cdragm,
1292	s dsens,
1293	s ycoefh,yu1,yv1)
1294	endif
1295
1296	CXXX Incrementation des flux
1297	DO i = 1, klon
1298	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
1299	fder(i) = dlw(i) + dsens(i)
1300	ENDDO
1301	CXXX
1302	IF (.not. turb_resolved) then !True only for LES
1303	DO k = 1, klev
1304	DO i = 1, klon
1305	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1306	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
1307	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1308	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
1309	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1310	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
1311	ENDDO
1312	ENDDO
1313	ENDIF
1314	C TRACEURS
1315
1316	if (iflag_trac.eq.1) then
1317	DO k = 1, klev
1318	DO i = 1, klon
1319	delp(i,k) = paprs(i,k)-paprs(i,k+1)
1320	ENDDO
1321	ENDDO
1322
1323	DO iq=1, nqmax
1324
1325	CALL cltrac(dtime,ycoefh,t_seri,
1326	s tr_seri(:,:,iq),source(:,iq),
1327	e paprs, pplay,delp,
1328	s d_tr_vdf(:,:,iq))
1329
1330	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
1331	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
1332
1333	ENDDO !nqmax
1334
1335	endif
1336
1337	IF (if_ebil.ge.2) THEN
1338	ztit='after clmain'
1339	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
1340	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1341	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1342	call diagphy(cell_area,ztit,ip_ebil
1343	e , zero_v, zero_v, zero_v, zero_v, sens
1344	e , zero_v, zero_v, zero_v, ztsol
1345	e , d_h_vcol, d_qt, d_ec
1346	s , fs_bound, fq_bound )
1347	END IF
1348	C
1349	c
1350	c Incrementer la temperature du sol
1351	c
1352	DO i = 1, klon
1353	ftsol(i) = ftsol(i) + d_ts(i)
1354	ENDDO
1355
1356	c Calculer la derive du flux infrarouge
1357	c
1358	DO i = 1, klon
1359	dlw(i) = - 4.0RSIGMAftsol(i)**3
1360	ENDDO
1361
1362	c-------------------------------
1363	endif ! fin du VENUS TEST
1364
1365	! tests: output tendencies
1366	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
1367	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
1368	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
1369	! call writefield_phy('physiq_d_ts',d_ts,1)
1370
1371	c===================================================================
1372	c Convection seche
1373	c===================================================================
1374	c
1375	d_t_ajs(:,:)=0.
1376	d_u_ajs(:,:)=0.
1377	d_v_ajs(:,:)=0.
1378	d_tr_ajs(:,:,:)=0.
1379	c
1380	IF(prt_level>9)WRITE(lunout,*)
1381	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
1382	s ,iflag_ajs
1383
1384	if(iflag_ajs.eq.0) then
1385	c Rien
1386	c ====
1387	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1388
1389	else if(iflag_ajs.eq.1) then
1390
1391	c Ajustement sec
1392	c ==============
1393	IF(prt_level>9)WRITE(lunout,*)'ajsec'
1394
1395	! ADAPTATION GCM POUR CP(T)
1396	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
1397	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
1398
1399	! ADAPTATION GCM POUR CP(T)
1400	do i=1,klon
1401	flux_ajs(i,1) = 0.0
1402	do j=2,klev
1403	flux_ajs(i,j) = flux_ajs(i,j-1)
1404	. + cpnew(i,j-1)/RG*d_t_ajs(i,j-1)/dtime
1405	. *(paprs(i,j-1)-paprs(i,j))
1406	enddo
1407	enddo
1408
1409	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1410	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1411	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1412	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1413	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1414	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1415
1416	if (iflag_trac.eq.1) then
1417	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1418	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1419	endif
1420	endif
1421
1422	! tests: output tendencies
1423	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
1424	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
1425	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
1426	c
1427	IF (if_ebil.ge.2) THEN
1428	ztit='after dry_adjust'
1429	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1430	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1431	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1432	call diagphy(cell_area,ztit,ip_ebil
1433	e , zero_v, zero_v, zero_v, zero_v, sens
1434	e , zero_v, zero_v, zero_v, ztsol
1435	e , d_h_vcol, d_qt, d_ec
1436	s , fs_bound, fq_bound )
1437	END IF
1438
1439	c====================================================================
1440	c RAYONNEMENT
1441	c====================================================================
1442	if (mod(itap,radpas) == 0) then
1443
1444	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1445
1446	! update mmean
1447	if (ok_chem) then
1448	mmean(:,:) = 0.
1449	do iq = 1,nqmax - nmicro
1450	mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq)
1451	enddo
1452	mmean(:,:) = 1./mmean(:,:)
1453	rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K
1454	endif
1455
1456	cc---------------------------------------------
1457	if (callnlte .or. callthermos) then
1458	if (ok_chem) then
1459
1460	! nlte : use computed chemical species
1461
1462	co2vmr_gcm(:,:) = tr_seri(:,:,i_co2)*mmean(:,:)/m_tr(i_co2)
1463	covmr_gcm(:,:) = tr_seri(:,:,i_co)*mmean(:,:)/m_tr(i_co)
1464	ovmr_gcm(:,:) = tr_seri(:,:,i_o)*mmean(:,:)/m_tr(i_o)
1465	n2vmr_gcm(:,:) = tr_seri(:,:,i_n2)*mmean(:,:)/m_tr(i_n2)
1466
1467	else
1468
1469	! nlte : use hedin climatology
1470
1471	call compo_hedin83_mod(pplay,rmu0,
1472	$ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1473	end if
1474	end if
1475
1476	c NLTE cooling from CO2 emission
1477	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1478
1479	IF(callnlte) THEN
1480	if(nltemodel.eq.0.or.nltemodel.eq.1) then
1481	c nltecool call not correct...
1482	c CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
1483	c $ tr_seri, d_t_nlte)
1484	abort_message='nltemodel=0 or 1 should not be used...'
1485	call abort_physic(modname,abort_message,1)
1486	else if(nltemodel.eq.2) then
1487	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1488	!! HEDIN instead of compo for this calculation
1489	! call compo_hedin83_mod(pplay,rmu0,
1490	! $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1491	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1492	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
1493	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1494	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
1495	if(ierr_nlte.gt.0) then
1496	write(,)
1497	$ 'WARNING: nlte_tcool output with error message',
1498	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
1499	write(,)'I will continue anyway'
1500	endif
1501	endif
1502
1503	ELSE
1504
1505	d_t_nlte(:,:)=0.
1506
1507	ENDIF
1508
1509	c Find number of layers for LTE radiation calculations
1510
1511	IF(callnlte .or. callnirco2)
1512	$ CALL nlthermeq(klon, klev, paprs, pplay)
1513
1514	cc---------------------------------------------
1515	c LTE radiative transfert / solar / IR matrix
1516	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1517	if (physideal) then
1518	CALL radlwsw_newtoncool(presnivs,t_seri)
1519	else
1520	CALL radlwsw
1521	e (dist, rmu0, fract, zzlev,
1522	e paprs, pplay,ftsol, t_seri)
1523	endif
1524
1525	c ALBEDO VARIATIONS: test for Yeon Joo Lee
1526	c increment to increase it for 20 Vd => +80%
1527	c heat(:,:)=heat(:,:)(1.+0.80((rjourvrai-356)+gmtime)/20.)
1528	c or to decrease it for 20 Vd => 1/1.8
1529	c heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.)
1530
1531	c ------------ ALBEDO VARIATIONS: scenarios for VCD
1532	c shape of relative variation from Lee et al 2019 (Fig 13b)
1533	c between 57 km (4e4 Pa) and 72 km (2.5e3 Pa), peak at 67 km (6e3 Pa)
1534	c do j=1,klev
1535	c factAlb = 0.
1536	c if ((presnivs(j).gt.6.e3).and.(presnivs(j).lt.4.e4)) then
1537	c factAlb = (log(presnivs(j))-log(4.e4))/(log(6.e3)-log(4.e4))
1538	c elseif ((presnivs(j).lt.6.e3).and.(presnivs(j).gt.2.5e3)) then
1539	c factAlb = (log(presnivs(j))-log(2.5e3))/(log(6.e3)-log(2.5e3))
1540	c endif
1541	c Increase by 50% (Minimum albedo)
1542	c heat(:,j)=heat(:,j)(1+factAlb0.5)
1543	c Decrease by 30% (Maximum albedo)
1544	c heat(:,j)=heat(:,j)(1-factAlb0.3)
1545	c enddo
1546	c ------------ END ALBEDO VARIATIONS
1547
1548	cc---------------------------------------------
1549	c CO2 near infrared absorption
1550	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1551
1552	d_t_nirco2(:,:)=0.
1553	if (callnirco2) then
1554	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1555	!! HEDIN instead of compo for this calculation
1556	! call compo_hedin83_mod(pplay,rmu0,
1557	! $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1558	! tr_hedin(:,:,:)=tr_seri(:,:,:)
1559	! tr_hedin(:,:,i_co2)=co2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_co2)
1560	! tr_hedin(:,:,i_co) = covmr_gcm(:,:)/mmean(:,:)*m_tr(i_co)
1561	! tr_hedin(:,:,i_o) = ovmr_gcm(:,:)/mmean(:,:)*m_tr(i_o)
1562	! tr_hedin(:,:,i_n2) = n2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_n2)
1563	! call nirco2abs (klon, klev, pplay, dist, nqmax, tr_hedin,
1564	! . rmu0, fract, d_t_nirco2)
1565	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1566	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
1567	. rmu0, fract, d_t_nirco2)
1568	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1569	endif
1570
1571
1572	cc---------------------------------------------
1573	c Net atmospheric radiative heating rate (K.s-1)
1574	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1575
1576	IF(callnlte.or.callnirco2) THEN
1577	CALL blendrad(klon, klev, pplay,heat,
1578	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
1579	ELSE
1580	dtsw(:,:)=heat(:,:)
1581	dtlw(:,:)=-1*cool(:,:)
1582	ENDIF
1583
1584	DO k=1,klev
1585	DO i=1,klon
1586	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
1587	ENDDO
1588	ENDDO
1589
1590
1591	cc---------------------------------------------
1592	c EUV heating rate (K.s-1)
1593	c ~~~~~~~~~~~~~~~~~~~~~~~~
1594
1595	d_t_euv(:,:)=0.
1596
1597	IF (callthermos) THEN
1598
1599	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
1600	$ rmu0,dtimerad,gmtime,
1601	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1602	!! HEDIN instead of compo for this calculation
1603	!! cf nlte_tcool and nirco2abs above !!
1604	! $ tr_hedin, d_tr, d_t_euv )
1605	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1606	$ tr_seri, d_tr, d_t_euv )
1607	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1608
1609	DO k=1,klev
1610	DO ig=1,klon
1611	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
1612	ENDDO
1613	ENDDO
1614
1615	ENDIF ! callthermos
1616
1617	ENDIF ! radpas
1618	c====================================================================
1619	c
1620	c Ajouter la tendance des rayonnements (tous les pas)
1621	c
1622	DO k = 1, klev
1623	DO i = 1, klon
1624	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
1625	ENDDO
1626	ENDDO
1627
1628	! increment physics counter
1629
1630	itap = itap + 1
1631	c====================================================================
1632
1633
1634	c==============================
1635	! -- MOLECULAR DIFFUSION ---
1636	c==============================
1637
1638	d_q_moldif(:,:,:)=0
1639
1640	IF (callthermos .and. ok_chem) THEN
1641
1642	call moldiff_red(klon, klev, nqmax,
1643	& pplay,paprs,t_seri, tr_seri, pdtphys,
1644	& d_t_euv,d_t_conduc,d_q_moldif)
1645
1646
1647	! --- update tendencies tracers ---
1648
1649	DO iq = 1, nqmax
1650	DO k=1,klev
1651	DO ig=1,klon
1652	tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
1653	& d_q_moldif(ig,k,iq)*dtime,1.e-30)
1654	ENDDO
1655	ENDDO
1656	ENDDO
1657
1658	ENDIF ! callthermos & ok_chem
1659
1660	c====================================================================
1661	! tests: output tendencies
1662	! call writefield_phy('physiq_dtrad',dtrad,klev)
1663
1664	IF (if_ebil.ge.2) THEN
1665	ztit='after rad'
1666	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1667	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1668	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1669	call diagphy(cell_area,ztit,ip_ebil
1670	e , topsw, toplw, solsw, sollw, zero_v
1671	e , zero_v, zero_v, zero_v, ztsol
1672	e , d_h_vcol, d_qt, d_ec
1673	s , fs_bound, fq_bound )
1674	END IF
1675	c
1676
1677	c====================================================================
1678	c Calcul des gravity waves FLOTT
1679	c====================================================================
1680	c
1681	c if (ok_orodr.or.ok_gw_nonoro) then
1682
1683	c CALCUL DE N2
1684	c UTILISE LA RELATION ENTRE N2 ET STABILITE
1685	c N2 = RG/T (dT/dz+RG/cp(T))
1686	c ET DONC EN N'UTILISE QUE LA TEMPERATURE, PAS teta.
1687
1688	do i=1,klon
1689	do k=2,klev
1690	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1691	enddo
1692	enddo
1693	do i=1,klon
1694	do k=2,klev
1695	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1696	zdtlev(i,k) = t_seri(i,k)-t_seri(i,k-1)
1697	zdzlev(i,k) = (zzlay(i,k)-zzlay(i,k-1))
1698	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
1699	. + RG/cpnew(i,k) )
1700	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1701	enddo
1702	zn2(i,1) = 1.e-12 ! securite
1703	enddo
1704
1705	c endif
1706
1707	c ----------------------------ORODRAG
1708	IF (ok_orodr) THEN
1709	c
1710	c selection des points pour lesquels le shema est actif:
1711	igwd=0
1712	DO i=1,klon
1713	itest(i)=0
1714	c IF ((zstd(i).gt.10.0)) THEN
1715	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1716	itest(i)=1
1717	igwd=igwd+1
1718	idx(igwd)=i
1719	ENDIF
1720	ENDDO
1721	c igwdim=MAX(1,igwd)
1722	c
1723	c A ADAPTER POUR VENUS!!! [ TN: c'est fait ! ]
1724	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
1725	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1726	e igwd,idx,itest,
1727	e t_seri, u_seri, v_seri,
1728	s zulow, zvlow, zustrdr, zvstrdr,
1729	s d_t_oro, d_u_oro, d_v_oro,
1730	s zublstrdr,zvblstrdr,znlow,zeff,zbl,
1731	s ztau,tau0,knu2,kbreak)
1732
1733	c print*,"d_u_oro=",d_u_oro(klon/2,:)
1734	c ajout des tendances
1735	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1736	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1737	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1738	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1739	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1740	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1741	c
1742	ELSE
1743	d_t_oro = 0.
1744	d_u_oro = 0.
1745	d_v_oro = 0.
1746	zustrdr = 0.
1747	zvstrdr = 0.
1748	zublstrdr = 0.
1749	zvblstrdr = 0.
1750	znlow = 0.
1751	zeff = 0.
1752	zbl = 0
1753	knu2 = 0
1754	kbreak = 0
1755	ztau = 0
1756	tau0 = 0.
1757	c
1758	ENDIF ! fin de test sur ok_orodr
1759	c
1760	! tests: output tendencies
1761	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
1762	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
1763	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
1764
1765	c ----------------------------OROLIFT
1766	IF (ok_orolf) THEN
1767	print*,"ok_orolf NOT IMPLEMENTED !"
1768	stop
1769	c
1770	c selection des points pour lesquels le shema est actif:
1771	igwd=0
1772	DO i=1,klon
1773	itest(i)=0
1774	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1775	itest(i)=1
1776	igwd=igwd+1
1777	idx(igwd)=i
1778	ENDIF
1779	ENDDO
1780	c igwdim=MAX(1,igwd)
1781	c
1782	c A ADAPTER POUR VENUS!!!
1783	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
1784	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1785	c e igwd,idx,itest,
1786	c e t_seri, u_seri, v_seri,
1787	c s zulow, zvlow, zustrli, zvstrli,
1788	c s d_t_lif, d_u_lif, d_v_lif )
1789
1790	c
1791	c ajout des tendances
1792	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1793	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1794	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1795	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1796	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1797	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1798	c
1799	ELSE
1800	d_t_lif = 0.
1801	d_u_lif = 0.
1802	d_v_lif = 0.
1803	zustrli = 0.
1804	zvstrli = 0.
1805	c
1806	ENDIF ! fin de test sur ok_orolf
1807
1808	c ---------------------------- NON-ORO GRAVITY WAVES
1809	IF(ok_gw_nonoro) then
1810
1811	! Obsolete
1812	! but used for VCD 1.1
1813	! call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1814	! e t_seri, u_seri, v_seri, paprs(klon/2+1,:),
1815	! o zustrhi,zvstrhi,
1816	! o d_t_hin, d_u_hin, d_v_hin)
1817
1818	! New routine based on Generic
1819	! used after VCD 1.1, for VCD 2.0
1820	call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs,
1821	e t_seri, u_seri, v_seri,
1822	o zustrhi,zvstrhi,
1823	o d_t_hin, d_u_hin, d_v_hin)
1824
1825	c ajout des tendances
1826
1827	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1828	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1829	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1830	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1831	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1832	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1833
1834	ELSE
1835	d_t_hin = 0.
1836	d_u_hin = 0.
1837	d_v_hin = 0.
1838	zustrhi = 0.
1839	zvstrhi = 0.
1840
1841	ENDIF ! fin de test sur ok_gw_nonoro
1842
1843	! tests: output tendencies
1844	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
1845	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
1846	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
1847
1848	c====================================================================
1849	c Transport de ballons
1850	c====================================================================
1851	if (ballons.eq.1) then
1852	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
1853	& latitude_deg,longitude_deg,
1854	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1855	C t,pplay,u,v,zphi) ! alt above planet average radius
1856	endif !ballons
1857
1858	c====================================================================
1859	c Bilan de mmt angulaire
1860	c====================================================================
1861	if (bilansmc.eq.1) then
1862	CMODDEB FLOTT
1863	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1864	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1865
1866	DO i = 1, klon
1867	zustrph(i)=0.
1868	zvstrph(i)=0.
1869	zustrcl(i)=0.
1870	zvstrcl(i)=0.
1871	ENDDO
1872	DO k = 1, klev
1873	DO i = 1, klon
1874	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1875	c (paprs(i,k)-paprs(i,k+1))/rg
1876	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1877	c (paprs(i,k)-paprs(i,k+1))/rg
1878	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1879	c (paprs(i,k)-paprs(i,k+1))/rg
1880	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1881	c (paprs(i,k)-paprs(i,k+1))/rg
1882	ENDDO
1883	ENDDO
1884
1885	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
1886	C ra,rg,romega,
1887	C latitude_deg,longitude_deg,pphis,
1888	C zustrdr,zustrli,zustrcl,
1889	C zvstrdr,zvstrli,zvstrcl,
1890	C paprs,u,v)
1891
1892	CCMODFIN FLOTT
1893	endif !bilansmc
1894
1895	c====================================================================
1896	c====================================================================
1897	c Calculer le transport de l'eau et de l'energie (diagnostique)
1898	c
1899	c A REVOIR POUR VENUS...
1900	c
1901	c CALL transp (paprs,ftsol,
1902	c e t_seri, q_seri, u_seri, v_seri, zphi,
1903	c s ve, vq, ue, uq)
1904	c
1905	c====================================================================
1906	c+jld ec_conser
1907	DO k = 1, klev
1908	DO i = 1, klon
1909	d_t_ec(i,k)=0.5/cpnew(i,k)
1910	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1911	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1912	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1913	END DO
1914	END DO
1915	do i=1,klon
1916	flux_ec(i,1) = 0.0
1917	do j=2,klev
1918	flux_ec(i,j) = flux_ec(i,j-1)
1919	. +cpnew(i,j-1)/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
1920	enddo
1921	enddo
1922
1923	c-jld ec_conser
1924	c====================================================================
1925	IF (if_ebil.ge.1) THEN
1926	ztit='after physic'
1927	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
1928	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1929	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1930	C Comme les tendances de la physique sont ajoute dans la dynamique,
1931	C on devrait avoir que la variation d'entalpie par la dynamique
1932	C est egale a la variation de la physique au pas de temps precedent.
1933	C Donc la somme de ces 2 variations devrait etre nulle.
1934	call diagphy(cell_area,ztit,ip_ebil
1935	e , topsw, toplw, solsw, sollw, sens
1936	e , zero_v, zero_v, zero_v, ztsol
1937	e , d_h_vcol, d_qt, d_ec
1938	s , fs_bound, fq_bound )
1939	C
1940	d_h_vcol_phy=d_h_vcol
1941	C
1942	END IF
1943	C
1944	c=======================================================================
1945	c SORTIES
1946	c=======================================================================
1947
1948	c Convertir les incrementations en tendances
1949	c
1950	DO k = 1, klev
1951	DO i = 1, klon
1952	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1953	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1954	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1955	ENDDO
1956	ENDDO
1957	c
1958	DO iq = 1, nqmax
1959	DO k = 1, klev
1960	DO i = 1, klon
1961	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1962	ENDDO
1963	ENDDO
1964	ENDDO
1965
1966	c mise à jour rho,mmean pour sorties
1967	if(callthermos) then
1968	call concentrations2(pplay,t_seri,tr_seri, nqmax)
1969	endif
1970
1971	c calcul vitesse verticale en m/s
1972	DO k = 1, klev
1973	DO i = 1, klon
1974	vertwind(i,k) = -omega(i,k)/(rho(i,k)*RG)
1975	END DO
1976	END DO
1977
1978	c------------------------
1979	c Calcul moment cinetique
1980	c------------------------
1981	c TEST VENUS...
1982	c mangtot = 0.0
1983	c DO k = 1, klev
1984	c DO i = 1, klon
1985	c mang(i,k) = RA*cos(latitude(i))
1986	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
1987	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
1988	c mangtot=mangtot+mang(i,k)
1989	c ENDDO
1990	c ENDDO
1991	c print*,"Moment cinetique total = ",mangtot
1992	c
1993	c------------------------
1994	c
1995	c Sauvegarder les valeurs de t et u a la fin de la physique:
1996	c
1997	DO k = 1, klev
1998	DO i = 1, klon
1999	u_ancien(i,k) = u_seri(i,k)
2000	t_ancien(i,k) = t_seri(i,k)
2001	ENDDO
2002	ENDDO
2003	c
2004	c=============================================================
2005	c Ecriture des sorties
2006	c=============================================================
2007	#ifndef MESOSCALE
2008	#ifdef CPP_IOIPSL
2009
2010	#ifdef histhf
2011	#include "write_histhf.h"
2012	#endif
2013
2014	#ifdef histday
2015	#include "write_histday.h"
2016	#endif
2017
2018	#ifdef histmth
2019	#include "write_histmth.h"
2020	#endif
2021
2022	#ifdef histins
2023	#include "write_histins.h"
2024	#endif
2025
2026	#endif
2027
2028	! XIOS outputs
2029	! This can be done ANYWHERE in the physics routines !
2030
2031	#ifdef CPP_XIOS
2032	! Send fields to XIOS: (NB these fields must also be defined as
2033	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
2034
2035	! 2D fields
2036
2037	CALL send_xios_field("phis",pphis)
2038	cell_area_out(:)=cell_area(:)
2039	if (is_north_pole_phy) cell_area_out(1)=cell_area(1)/nbp_lon
2040	if (is_south_pole_phy) cell_area_out(klon)=cell_area(klon)/nbp_lon
2041	CALL send_xios_field("aire",cell_area_out)
2042	CALL send_xios_field("tsol",ftsol)
2043	CALL send_xios_field("psol",paprs(:,1))
2044	CALL send_xios_field("cdragh",cdragh)
2045	CALL send_xios_field("cdragm",cdragm)
2046
2047	CALL send_xios_field("tops",topsw)
2048	CALL send_xios_field("topl",toplw)
2049	CALL send_xios_field("sols",solsw)
2050	CALL send_xios_field("soll",sollw)
2051
2052	! 3D fields
2053
2054	CALL send_xios_field("temp",t_seri)
2055	CALL send_xios_field("pres",pplay)
2056	CALL send_xios_field("geop",zphi)
2057	CALL send_xios_field("vitu",u_seri)
2058	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2059	CALL send_xios_field("vitv",-1.*v_seri)
2060	CALL send_xios_field("vitw",omega)
2061	CALL send_xios_field("vitwz",vertwind)
2062	CALL send_xios_field("Kz",ycoefh)
2063	CALL send_xios_field("mmean",mmean)
2064	CALL send_xios_field("rho",rho)
2065	CALL send_xios_field("BV2",zn2)
2066
2067	CALL send_xios_field("dudyn",d_u_dyn)
2068	CALL send_xios_field("duvdf",d_u_vdf)
2069	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2070	CALL send_xios_field("dvvdf",-1.*d_v_vdf)
2071	CALL send_xios_field("duajs",d_u_ajs)
2072	CALL send_xios_field("dugwo",d_u_oro)
2073	CALL send_xios_field("dugwno",d_u_hin)
2074	CALL send_xios_field("dvgwno",-1.*d_v_hin)
2075	CALL send_xios_field("dumolvis",d_u_molvis)
2076	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
2077	CALL send_xios_field("dvmolvis",-1.*d_v_molvis)
2078	CALL send_xios_field("dtdyn",d_t_dyn)
2079	CALL send_xios_field("dtphy",d_t)
2080	CALL send_xios_field("dtvdf",d_t_vdf)
2081	CALL send_xios_field("dtajs",d_t_ajs)
2082	CALL send_xios_field("dtswr",dtsw)
2083	CALL send_xios_field("dtswrNLTE",d_t_nirco2)
2084	CALL send_xios_field("dtswrLTE",heat)
2085	CALL send_xios_field("dtlwr",dtlw)
2086	CALL send_xios_field("dtlwrNLTE",d_t_nlte)
2087	CALL send_xios_field("dtlwrLTE",-1.*cool)
2088	CALL send_xios_field("dteuv",d_t_euv)
2089	CALL send_xios_field("dtcond",d_t_conduc)
2090	CALL send_xios_field("dtec",d_t_ec)
2091
2092	CALL send_xios_field("SWnet",swnet(:,1:klev))
2093	CALL send_xios_field("LWnet",lwnet(:,1:klev))
2094	CALL send_xios_field("fluxvdf",fluxt)
2095	CALL send_xios_field("fluxdyn",flux_dyn)
2096	CALL send_xios_field("fluxajs",flux_ajs)
2097	CALL send_xios_field("fluxec",flux_ec)
2098
2099	! when using tracers
2100
2101	if (iflag_trac == 1) then
2102
2103	! production and destruction rate, cm-3.s-1
2104	! Beware of the context*.xml file !!
2105	if ((tr_scheme == 3) .and. (ok_chem)) THEN
2106	do iq = 1,nqmax - nmicro
2107	if ((iq.eq.i_o).or.(iq.eq.i_co)) THEN
2108	call send_xios_field("prod_"//tname(iq),
2109	$ prod_tr(:,:,iq))
2110	call send_xios_field("loss_"//tname(iq),
2111	$ loss_tr(:,:,iq))
2112	end if
2113	!if (iq.eq.i_o) then
2114	! call send_xios_field('prod_o', prod_tr(:,:,iq))
2115	! call send_xios_field('loss_o', loss_tr(:,:,iq))
2116	!end if
2117	!if (iq.eq.i_co) then
2118	! call send_xios_field('prod_co', prod_tr(:,:,iq))
2119	! call send_xios_field('loss_co', loss_tr(:,:,iq))
2120	!end if
2121	end do
2122	end if
2123
2124	! tracers in gas phase, volume mixing ratio
2125
2126	do iq = 1,nqmax - nmicro
2127	call send_xios_field(tname(iq),
2128	$ tr_seri(:,:,iq)*mmean(:,:)/m_tr(iq))
2129	end do
2130
2131	! tracers in gas phase, column densities
2132
2133	do iq = 1,nqmax - nmicro
2134	col_dens_tr(:,iq)=0.
2135	if (type_tr(iq).eq.1) THEN
2136	do k = 1, klev
2137	col_dens_tr(:,iq) = col_dens_tr(:,iq) +
2138	$ tr_seri(:,k,iq)* (paprs(:,k)-paprs(:,k+1)) / RG
2139	end do
2140	call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq))
2141	end if
2142	end do
2143
2144	! tracers in liquid phase, volume mixing ratio
2145
2146	if ((tr_scheme == 3) .and. (cl_scheme == 1)) THEN
2147	call send_xios_field(tname(i_h2oliq),
2148	$ tr_seri(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
2149	call send_xios_field(tname(i_h2so4liq),
2150	$ tr_seri(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
2151	if (ok_sedim) then
2152	call send_xios_field("Fsedim",fsedim(:,1:klev))
2153	end if
2154	end if
2155
2156	! aeronomical emissions
2157
2158	call send_xios_field("no_emis",no_emission)
2159	call send_xios_field("o2_emis",o2_emission)
2160
2161	! chemical iterations
2162
2163	if (tr_scheme.eq.3) call send_xios_field("iter",real(iter))
2164
2165	end if
2166
2167	IF (callthermos .and. ok_chem) THEN
2168	CALL send_xios_field("d_qmoldifCO2",d_q_moldif(:,:,i_co2))
2169	CALL send_xios_field("d_qmoldifO3p",d_q_moldif(:,:,i_o))
2170	CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
2171	ENDIF
2172
2173	!! DEBUG
2174	! if (is_master) print*,"itauphy=",itap
2175	! if (itap.eq.10) lafin=.true.
2176
2177	if (lafin.and.is_omp_master) then
2178	write(,) "physiq: call xios_context_finalize"
2179	call xios_context_finalize
2180	endif
2181
2182	#endif
2183	#else
2184	! Outputs MESOSCALE
2185	CALL allocate_comm_wrf(klon,klev)
2186	comm_HR_SW(1:klon,1:klev) = dtsw(1:klon,1:klev)
2187	comm_HR_LW(1:klon,1:klev) = dtlw(1:klon,1:klev)
2188	comm_DT_RAD(1:klon,1:klev) = d_t_rad(1:klon,1:klev)
2189	IF (turb_resolved) THEN
2190	open(17,file='hrdyn.txt',form='formatted',status='old')
2191	rewind(17)
2192	DO k=1,klev
2193	read(17,*) dt_dyn(k)
2194	ENDDO
2195	close(17)
2196
2197	do i=1,klon
2198	d_t(i,:)=d_t(i,:)+dt_dyn(:)
2199	comm_HR_DYN(i,:) = dt_dyn(:)
2200	enddo
2201	ELSE
2202	comm_HR_DYN(1:klon,1:klev) = d_t_dyn(1:klon,1:klev)
2203	comm_DT_VDF(1:klon,1:klev) = d_t_vdf(1:klon,1:klev)
2204	comm_DT_AJS(1:klon,1:klev) = d_t_ajs(1:klon,1:klev)
2205	ENDIF
2206	comm_DT(1:klon,1:klev)=d_t(1:klon,1:klev)
2207	#endif
2208
2209
2210	c====================================================================
2211	c Si c'est la fin, il faut conserver l'etat de redemarrage
2212	c====================================================================
2213	c
2214	IF (lafin) THEN
2215	itau_phy = itau_phy + itap
2216	CALL phyredem ("restartphy.nc")
2217
2218	c--------------FLOTT
2219	CMODEB LOTT
2220	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
2221	C DU BILAN DE MOMENT ANGULAIRE.
2222	if (bilansmc.eq.1) then
2223	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
2224	close(27)
2225	close(28)
2226	endif !bilansmc
2227	CMODFIN
2228	c-------------
2229	c--------------SLEBONNOIS
2230	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
2231	C DES BALLONS
2232	if (ballons.eq.1) then
2233	write(,)'Fermeture des ballons*.out'
2234	close(30)
2235	close(31)
2236	close(32)
2237	close(33)
2238	close(34)
2239	endif !ballons
2240	c-------------
2241	ENDIF
2242
2243	END SUBROUTINE physiq
2244
2245	END MODULE physiq_mod
2246

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

Download in other formats:

Original Format