Context Navigation

physiq_mod.F @ 3900

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

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

Download in other formats:

Original Format