Context Navigation

physiq_mod.F @ 3656

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

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

Download in other formats:

Original Format