Context Navigation

physiq_mod.F @ 1704

Last change on this file since 1704 was 1691, checked in by slebonnois, 8 years ago
SL: corrections after testing of cloud microphysics in 1D
File size: 59.1 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, plevmoy, tmoy,
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 mod_phys_lmdz_para, only : is_parallel,jj_nb,
68	& is_north_pole_phy,
69	& is_south_pole_phy
70	USE phys_state_var_mod ! Variables sauvegardees de la physique
71	USE write_field_phy
72	USE iophy
73	USE cpdet_phy_mod, only: cpdet, t2tpot
74	USE chemparam_mod
75	USE conc
76	USE compo_hedin83_mod2
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 mod_phys_lmdz_omp_data, ONLY: is_omp_master
83	#ifdef CPP_XIOS
84	use xios_output_mod, only: initialize_xios_output,
85	& update_xios_timestep,
86	& send_xios_field
87	use wxios, only: wxios_context_init, xios_context_finalize
88	#endif
89	IMPLICIT none
90	c======================================================================
91	c CLEFS CPP POUR LES IO
92	c =====================
93	c#define histhf
94	#define histday
95	#define histmth
96	#define histins
97	c======================================================================
98	#include "dimsoil.h"
99	#include "clesphys.h"
100	#include "iniprint.h"
101	#include "timerad.h"
102	#include "tabcontrol.h"
103	#include "nirdata.h"
104	#include "hedin.h"
105	c======================================================================
106	LOGICAL ok_journe ! sortir le fichier journalier
107	save ok_journe
108	c PARAMETER (ok_journe=.true.)
109	c
110	LOGICAL ok_mensuel ! sortir le fichier mensuel
111	save ok_mensuel
112	c PARAMETER (ok_mensuel=.true.)
113	c
114	LOGICAL ok_instan ! sortir le fichier instantane
115	save ok_instan
116	c PARAMETER (ok_instan=.true.)
117	c
118	c======================================================================
119	c
120	c Variables argument:
121	c
122	INTEGER nlon
123	INTEGER nlev
124	INTEGER nqmax
125	REAL rjourvrai
126	REAL gmtime
127	REAL pdtphys
128	LOGICAL debut, lafin
129	REAL paprs(klon,klev+1)
130	REAL pplay(klon,klev)
131	REAL pphi(klon,klev)
132	REAL pphis(klon)
133	REAL presnivs(klev)
134
135	! ADAPTATION GCM POUR CP(T)
136	REAL ppk(klon,klev)
137
138	REAL u(klon,klev)
139	REAL v(klon,klev)
140	REAL t(klon,klev)
141	REAL qx(klon,klev,nqmax)
142
143	REAL d_u_dyn(klon,klev)
144	REAL d_t_dyn(klon,klev)
145
146	REAL flxmw(klon,klev)
147
148	! A VIRER POUR DYNAMICO
149	REAL,INTENT(IN) :: plevmoy(klev+1) ! planet-averaged mean pressure (Pa) at interfaces
150	REAL,INTENT(IN) :: tmoy(klev) ! planet-averaged mean temperature (Pa) at mid-layers
151
152	REAL d_u(klon,klev)
153	REAL d_v(klon,klev)
154	REAL d_t(klon,klev)
155	REAL d_qx(klon,klev,nqmax)
156	REAL d_ps(klon)
157
158	logical ok_hf
159	real ecrit_hf
160	integer nid_hf
161	save ok_hf, ecrit_hf, nid_hf
162
163	#ifdef histhf
164	data ok_hf,ecrit_hf/.true.,0.25/
165	#else
166	data ok_hf/.false./
167	#endif
168
169	c Variables propres a la physique
170	c
171	INTEGER,save :: itap ! compteur pour la physique
172	REAL delp(klon,klev) ! epaisseur d'une couche
173	REAL omega(klon,klev) ! vitesse verticale en Pa/s
174
175
176	INTEGER igwd,idx(klon),itest(klon)
177	c
178	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
179
180	REAL zulow(klon),zvlow(klon)
181	REAL zustrdr(klon), zvstrdr(klon)
182	REAL zustrli(klon), zvstrli(klon)
183	REAL zustrhi(klon), zvstrhi(klon)
184
185	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
186	real zdtlev(klon,klev),zdzlev(klon,klev)
187	real ztlev(klon,klev)
188	real zn2(klon,klev) ! BV^2 at plev
189
190	c Pour les bilans de moment angulaire,
191	integer bilansmc
192	c Pour le transport de ballons
193	integer ballons
194	c j'ai aussi besoin
195	c du stress de couche limite a la surface:
196
197	REAL zustrcl(klon),zvstrcl(klon)
198
199	c et du stress total c de la physique:
200
201	REAL zustrph(klon),zvstrph(klon)
202
203	c Variables locales:
204	c
205	REAL cdragh(klon) ! drag coefficient pour T and Q
206	REAL cdragm(klon) ! drag coefficient pour vent
207	c
208	cAA Pour TRACEURS
209	cAA
210	REAL,save,allocatable :: source(:,:)
211	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
212	REAL yu1(klon) ! vents dans la premiere couche U
213	REAL yv1(klon) ! vents dans la premiere couche V
214
215	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
216	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
217	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
218	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
219	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
220	c
221	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
222
223	c======================================================================
224	c
225	c Declaration des procedures appelees
226	c
227	EXTERNAL ajsec ! ajustement sec
228	EXTERNAL clmain ! couche limite
229	EXTERNAL hgardfou ! verifier les temperatures
230	c EXTERNAL orbite ! calculer l'orbite
231	EXTERNAL phyetat0 ! lire l'etat initial de la physique
232	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
233	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
234	EXTERNAL suphec ! initialiser certaines constantes
235	EXTERNAL transp ! transport total de l'eau et de l'energie
236	EXTERNAL abort_gcm
237	EXTERNAL printflag
238	EXTERNAL zenang
239	EXTERNAL diagetpq
240	EXTERNAL conf_phys
241	EXTERNAL diagphy
242	EXTERNAL mucorr
243	EXTERNAL nirco2abs
244	EXTERNAL nir_leedat
245	EXTERNAL nltecool
246	EXTERNAL nlte_tcool
247	EXTERNAL nlte_setup
248	EXTERNAL blendrad
249	EXTERNAL nlthermeq
250	EXTERNAL euvheat
251	EXTERNAL param_read
252	EXTERNAL param_read_e107
253	EXTERNAL conduction
254	EXTERNAL molvis
255	EXTERNAL moldiff_red
256
257	c
258	c Variables locales
259	c
260	CXXX PB
261	REAL fluxt(klon,klev) ! flux turbulent de chaleur
262	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
263	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
264	c
265	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
266	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
267	REAL flux_ec(klon,klev) ! flux de chaleur Ec
268	c
269	REAL tmpout(klon,klev) ! K s-1
270
271	INTEGER itaprad
272	SAVE itaprad
273	c
274	REAL dist, rmu0(klon), fract(klon)
275	REAL zdtime, zlongi
276	c
277	INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev
278	c
279	REAL zphi(klon,klev)
280	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
281	real tsurf(klon)
282
283	c va avec nlte_tcool
284	INTEGER ierr_nlte
285	REAL varerr
286
287	c Variables du changement
288	c
289	c ajs: ajustement sec
290	c vdf: couche limite (Vertical DiFfusion)
291	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
292	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
293	c
294	REAL d_ts(klon)
295	c
296	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
297	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
298	c
299	CMOD LOTT: Tendances Orography Sous-maille
300	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
301	REAL d_t_oro(klon,klev)
302	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
303	REAL d_t_lif(klon,klev)
304	C Tendances Ondes de G non oro (runs strato).
305	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
306	REAL d_t_hin(klon,klev)
307
308	c Tendencies due to radiative scheme [K/s]
309	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
310	c are not computed at each physical timestep
311	c therefore, they are defined and saved in phys_state_var_mod
312
313	c Tendencies due to molecular viscosity and conduction
314	real d_t_conduc(klon,klev) ! [K/s]
315	real d_u_molvis(klon,klev) ! (m/s) /s
316	real d_v_molvis(klon,klev) ! (m/s) /s
317
318	c Tendencies due to molecular diffusion
319	real d_q_moldif(klon,klev,nqmax)
320
321	c
322	c Variables liees a l'ecriture de la bande histoire physique
323	c
324	INTEGER ecrit_mth
325	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
326	c
327	INTEGER ecrit_day
328	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
329	c
330	INTEGER ecrit_ins
331	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
332	c
333	integer itau_w ! pas de temps ecriture = itap + itau_phy
334
335	c Variables locales pour effectuer les appels en serie
336	c
337	REAL t_seri(klon,klev)
338	REAL u_seri(klon,klev), v_seri(klon,klev)
339	c
340	REAL :: tr_seri(klon,klev,nqmax)
341	REAL :: d_tr(klon,klev,nqmax)
342
343	c Variables tendance sedimentation
344
345	REAL :: m0_mode1(klev,2),m0_mode2(klev,2)
346	REAL :: m3_mode1(klev,3),m3_mode2 (klev,3)
347	REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2)
348	REAL :: aer_flux(klev)
349	REAL :: d_tr_sed(klon,klev,nqmax)
350	REAL :: d_tr_ssed(klon)
351	c
352	c pour ioipsl
353	INTEGER nid_day, nid_mth, nid_ins
354	SAVE nid_day, nid_mth, nid_ins
355	INTEGER nhori, nvert, idayref
356	REAL zsto, zout, zsto1, zsto2, zero
357	parameter (zero=0.0e0)
358	real zjulian
359	save zjulian
360
361	CHARACTER*2 str2
362	character*20 modname
363	character*80 abort_message
364	logical ok_sync
365
366	character*30 nom_fichier
367	character*10 varname
368	character*40 vartitle
369	character*20 varunits
370	C Variables liees au bilan d'energie et d'enthalpi
371	REAL ztsol(klon)
372	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
373	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
374	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
375	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
376	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
377	REAL d_h_vcol_phy
378	REAL fs_bound, fq_bound
379	SAVE d_h_vcol_phy
380	REAL zero_v(klon),zero_v2(klon,klev)
381	CHARACTER*15 ztit
382	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
383	SAVE ip_ebil
384	DATA ip_ebil/2/
385	INTEGER if_ebil ! level for energy conserv. dignostics
386	SAVE if_ebil
387	c+jld ec_conser
388	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
389	c-jld ec_conser
390
391	c TEST VENUS...
392	REAL mang(klon,klev) ! moment cinetique
393	REAL mangtot ! moment cinetique total
394
395	c cell_area for outputs in hist*
396	REAL cell_area_out(klon)
397
398	c Declaration des constantes et des fonctions thermodynamiques
399	c
400	#include "YOMCST.h"
401
402	c======================================================================
403	c INITIALISATIONS
404	c================
405
406	modname = 'physiq'
407	ok_sync=.TRUE.
408
409	bilansmc = 0
410	ballons = 0
411	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
412	if (is_parallel) then
413	bilansmc = 0
414	ballons = 0
415	endif
416
417	IF (if_ebil.ge.1) THEN
418	DO i=1,klon
419	zero_v(i)=0.
420	END DO
421	DO i=1,klon
422	DO j=1,klev
423	zero_v2(i,j)=0.
424	END DO
425	END DO
426	END IF
427
428	c PREMIER APPEL SEULEMENT
429	c========================
430	IF (debut) THEN
431	allocate(source(klon,nqmax))
432
433	#ifdef CPP_XIOS
434	! Initialize XIOS context
435	write(,) "physiq: call wxios_context_init"
436	CALL wxios_context_init
437	#endif
438
439	CALL suphec ! initialiser constantes et parametres phys.
440
441	IF (if_ebil.ge.1) d_h_vcol_phy=0.
442	c
443	c appel a la lecture du physiq.def
444	c
445	call conf_phys(ok_journe, ok_mensuel,
446	. ok_instan,
447	. if_ebil)
448
449	call phys_state_var_init
450	c
451	c Initialising Hedin model for upper atm
452	c (to be revised when coupled to chemistry) :
453	call conc_init
454	c
455	c Initialiser les compteurs:
456	c
457	itap = 0
458	itaprad = 0
459	c
460	c Lecture startphy.nc :
461	c
462	CALL phyetat0 ("startphy.nc")
463
464	c dtime est defini dans tabcontrol.h et lu dans startphy
465	c pdtphys est calcule a partir des nouvelles conditions:
466	c Reinitialisation du pas de temps physique quand changement
467	IF (ABS(dtime-pdtphys).GT.0.001) THEN
468	WRITE(lunout,*) 'Pas physique a change',dtime,
469	. pdtphys
470	c abort_message='Pas physique n est pas correct '
471	c call abort_gcm(modname,abort_message,1)
472	c----------------
473	c pour initialiser convenablement le time_counter, il faut tenir compte
474	c du changement de dtime en changeant itau_phy (point de depart)
475	itau_phy = NINT(itau_phy*dtime/pdtphys)
476	c----------------
477	dtime=pdtphys
478	ENDIF
479
480	radpas = NINT( RDAY/pdtphys/nbapp_rad)
481
482	CALL printflag( ok_journe,ok_instan )
483
484	#ifdef CPP_XIOS
485
486	write(,) "physiq: call initialize_xios_output"
487	call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY,
488	& presnivs,pseudoalt)
489	#endif
490
491	c
492	c---------
493	c FLOTT
494	IF (ok_orodr) THEN
495	DO i=1,klon
496	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
497	ENDDO
498	CALL SUGWD(klon,klev,paprs,pplay)
499	DO i=1,klon
500	zuthe(i)=0.
501	zvthe(i)=0.
502	if(zstd(i).gt.10.)then
503	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
504	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
505	endif
506	ENDDO
507	ENDIF
508
509	if (bilansmc.eq.1) then
510	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
511	C DU BILAN DE MOMENT ANGULAIRE.
512	open(27,file='aaam_bud.out',form='formatted')
513	open(28,file='fields_2d.out',form='formatted')
514	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
515	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
516	endif !bilansmc
517
518	c--------------SLEBONNOIS
519	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
520	C DES BALLONS
521	if (ballons.eq.1) then
522	open(30,file='ballons-lat.out',form='formatted')
523	open(31,file='ballons-lon.out',form='formatted')
524	open(32,file='ballons-u.out',form='formatted')
525	open(33,file='ballons-v.out',form='formatted')
526	open(34,file='ballons-alt.out',form='formatted')
527	write(,)'Ouverture des ballons*.out'
528	endif !ballons
529	c-------------
530
531	c---------
532	C TRACEURS
533	C source dans couche limite
534	source(:,:) = 0.0 ! pas de source, pour l'instant
535	c---------
536
537	c---------
538	c INITIALIZE THERMOSPHERIC PARAMETERS
539	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
540
541	if (callthermos) then
542	if(solvarmod.eq.0) call param_read
543	if(solvarmod.eq.1) call param_read_e107
544	endif
545
546	c Initialisation (recomputed in concentration2)
547	do ig=1,klon
548	do j=1,klev
549	rnew(ig,j)=R
550	cpnew(ig,j)=cpdet(t(ig,j))
551	mmean(ig,j)=RMD
552	akknew(ig,j)=1.e-4
553	rho(ig,j)=pplay(ig,j)mmean(ig,j)1e-3/(rnew(ig,j)*t(ig,j))
554	enddo
555	c stop
556
557	enddo
558
559	IF(callthermos.or.callnlte.or.callnirco2) THEN
560	call compo_hedin83_init2
561	ENDIF
562	if (callnlte.and.nltemodel.eq.2) call nlte_setup
563	if (callnirco2.and.nircorr.eq.1) call nir_leedat
564	c---------
565
566	c
567	c Verifications:
568	c
569	IF (nlon .NE. klon) THEN
570	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
571	. klon
572	abort_message='nlon et klon ne sont pas coherents'
573	call abort_gcm(modname,abort_message,1)
574	ENDIF
575	IF (nlev .NE. klev) THEN
576	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
577	. klev
578	abort_message='nlev et klev ne sont pas coherents'
579	call abort_gcm(modname,abort_message,1)
580	ENDIF
581	c
582	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
583	$ THEN
584	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
585	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
586	abort_message='Nbre d appels au rayonnement insuffisant'
587	call abort_gcm(modname,abort_message,1)
588	ENDIF
589	c
590	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
591	. iflag_ajs
592	c
593	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
594	IF (ok_mensuel) THEN
595	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
596	. ecrit_mth
597	ENDIF
598
599	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
600	IF (ok_journe) THEN
601	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
602	. ecrit_day
603	ENDIF
604
605	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
606	IF (ok_instan) THEN
607	WRITE(lunout,*)'La frequence de sortie instant. est de ',
608	. ecrit_ins
609	ENDIF
610
611	c Initialisation des sorties
612	c===========================
613
614	#ifdef CPP_IOIPSL
615
616	#ifdef histhf
617	#include "ini_histhf.h"
618	#endif
619
620	#ifdef histday
621	#include "ini_histday.h"
622	#endif
623
624	#ifdef histmth
625	#include "ini_histmth.h"
626	#endif
627
628	#ifdef histins
629	#include "ini_histins.h"
630	#endif
631
632	#endif
633
634	c
635	c Initialiser les valeurs de u pour calculs tendances
636	c (pour T, c'est fait dans phyetat0)
637	c
638	DO k = 1, klev
639	DO i = 1, klon
640	u_ancien(i,k) = u(i,k)
641	ENDDO
642	ENDDO
643
644	c---------
645	c Ecriture fichier initialisation
646	c PRINT*,'Ecriture Initial_State.csv'
647	c OPEN(88,file='Trac_Point.csv',
648	c & form='formatted')
649	c---------
650
651	c---------
652	c Initialisation des parametres des nuages
653	c===============================================
654
655	c MICROPHY SANS CHIMIE: seulement si full microphy (cl_scheme=2)
656
657	if (ok_chem.and..not.ok_cloud) then
658	print*,"LA CHIMIE A BESOIN DE LA MICROPHYSIQUE"
659	print*,"ok_cloud doit etre = a ok_chem"
660	stop
661	endif
662	if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.1)) then
663	print*,"cl_scheme=1 doesnot work without chemistry"
664	stop
665	endif
666	if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.2)) then
667	print*,"Full microphysics without chemistry"
668	c indexation of microphys tracers
669	CALL chemparam_ini()
670	endif
671
672	c number of microphysical tracers
673	nmicro=0
674	if (ok_cloud .AND. (cl_scheme.eq.1)) nmicro=2
675	if (ok_cloud .AND. (cl_scheme.eq.2)) nmicro=12
676
677	c CAS 1D POUR MICROPHYS Aurelien
678	if ((nlon .EQ. 1) .AND. ok_cloud .AND. (cl_scheme.eq.1)) then
679	PRINT*,'Open profile_cloud_parameters.csv'
680	OPEN(66,file='profile_cloud_parameters.csv',
681	& form='formatted')
682	endif
683
684	if ((nlon .EQ. 1) .AND. ok_sedim .AND. (cl_scheme.eq.1)) then
685	PRINT*,'Open profile_cloud_sedim.csv'
686	OPEN(77,file='profile_cloud_sedim.csv',
687	& form='formatted')
688	endif
689
690	c INIT PHOTOCHEMISTRY ! includes the indexation of microphys tracers
691	if ((nlon .GT. 1) .AND. ok_chem) then
692	c !!! DONC 3D !!!
693	CALL chemparam_ini()
694	endif
695
696	c INIT MICROPHYS SCHEME 1 (AURELIEN)
697	if ((nlon .GT. 1) .AND. ok_cloud .AND. (cl_scheme.eq.1)) then
698	c !!! DONC 3D !!!
699	CALL cloud_ini(nlon,nlev)
700	endif
701
702	ENDIF ! debut
703	c======================================================================
704	c======================================================================
705	! ------------------------------------------------------
706	! Initializations done at every physical timestep:
707	! ------------------------------------------------------
708
709	c Mettre a zero des variables de sortie (pour securite)
710	c
711	DO i = 1, klon
712	d_ps(i) = 0.0
713	ENDDO
714	DO k = 1, klev
715	DO i = 1, klon
716	d_t(i,k) = 0.0
717	d_u(i,k) = 0.0
718	d_v(i,k) = 0.0
719	ENDDO
720	ENDDO
721	DO iq = 1, nqmax
722	DO k = 1, klev
723	DO i = 1, klon
724	d_qx(i,k,iq) = 0.0
725	ENDDO
726	ENDDO
727	ENDDO
728	c
729	c Ne pas affecter les valeurs entrees de u, v, h, et q
730	c
731	DO k = 1, klev
732	DO i = 1, klon
733	t_seri(i,k) = t(i,k)
734	u_seri(i,k) = u(i,k)
735	v_seri(i,k) = v(i,k)
736	ENDDO
737	ENDDO
738	DO iq = 1, nqmax
739	DO k = 1, klev
740	DO i = 1, klon
741	tr_seri(i,k,iq) = qx(i,k,iq)
742	ENDDO
743	ENDDO
744	ENDDO
745	C
746	DO i = 1, klon
747	ztsol(i) = ftsol(i)
748	ENDDO
749	C
750	IF (if_ebil.ge.1) THEN
751	ztit='after dynamic'
752	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
753	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
754	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
755	C Comme les tendances de la physique sont ajoute dans la dynamique,
756	C on devrait avoir que la variation d'entalpie par la dynamique
757	C est egale a la variation de la physique au pas de temps precedent.
758	C Donc la somme de ces 2 variations devrait etre nulle.
759	call diagphy(cell_area,ztit,ip_ebil
760	e , zero_v, zero_v, zero_v, zero_v, zero_v
761	e , zero_v, zero_v, zero_v, ztsol
762	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
763	s , fs_bound, fq_bound )
764	END IF
765
766	c====================================================================
767	c XIOS outputs
768
769	#ifdef CPP_XIOS
770	! update XIOS time/calendar
771	call update_xios_timestep
772	#endif
773
774	c====================================================================
775	c Diagnostiquer la tendance dynamique
776	c
777	IF (ancien_ok) THEN
778	DO k = 1, klev
779	DO i = 1, klon
780	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
781	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
782	ENDDO
783	ENDDO
784
785	! ADAPTATION GCM POUR CP(T)
786	do i=1,klon
787	flux_dyn(i,1) = 0.0
788	do j=2,klev
789	flux_dyn(i,j) = flux_dyn(i,j-1)
790	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
791	enddo
792	enddo
793
794	ELSE
795	DO k = 1, klev
796	DO i = 1, klon
797	d_u_dyn(i,k) = 0.0
798	d_t_dyn(i,k) = 0.0
799	ENDDO
800	ENDDO
801	ancien_ok = .TRUE.
802	ENDIF
803	c====================================================================
804
805	c Calcule de vitesse verticale a partir de flux de masse verticale
806	DO k = 1, klev
807	DO i = 1, klon
808	omega(i,k) = RG*flxmw(i,k) / cell_area(i)
809	END DO
810	END DO
811
812	c
813	c Ajouter le geopotentiel du sol:
814	c
815	DO k = 1, klev
816	DO i = 1, klon
817	zphi(i,k) = pphi(i,k) + pphis(i)
818	ENDDO
819	ENDDO
820
821	c calcul du geopotentiel aux niveaux intercouches
822	c ponderation des altitudes au niveau des couches en dp/p
823
824	DO k=1,klev
825	DO i=1,klon
826	zzlay(i,k)=zphi(i,k)/RG ! [m]
827	ENDDO
828	ENDDO
829	DO i=1,klon
830	zzlev(i,1)=pphis(i)/RG ! [m]
831	ENDDO
832	DO k=2,klev
833	DO i=1,klon
834	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
835	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
836	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
837	ENDDO
838	ENDDO
839	DO i=1,klon
840	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
841	ENDDO
842
843	c====================================================================
844	c
845	c Verifier les temperatures
846	c
847	CALL hgardfou(t_seri,ftsol,'debutphy')
848	c====================================================================
849	c
850	c Incrementer le compteur de la physique
851	c
852	itap = itap + 1
853
854	c====================================================================
855	c Orbite et eclairement
856	c====================================================================
857
858	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
859	c donc pas de variations de Ls, ni de dist.
860	c La seule chose qui compte, c'est la rotation de la planete devant
861	c le Soleil...
862
863	zlongi = 0.0
864	dist = 0.72 ! en UA
865
866	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
867	c a sa valeur, et prendre en compte leur evolution,
868	c il faudra refaire un orbite.F...
869	c CALL orbite(zlongi,dist)
870
871	IF (cycle_diurne) THEN
872	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
873	CALL zenang(zlongi,gmtime,zdtime,latitude_deg,longitude_deg,
874	& rmu0,fract)
875	ELSE
876	call mucorr(klon,zlongi,latitude_deg,rmu0,fract)
877	ENDIF
878
879	c====================================================================
880	c Calcul des tendances traceurs
881	c====================================================================
882
883	if (iflag_trac.eq.1) then
884
885	if (tr_scheme.eq.1) then
886	! Case 1: pseudo-chemistry with relaxation toward fixed profile
887
888	call phytrac_relax (debut,lafin,nqmax,
889	I nlon,nlev,dtime,pplay,
890	O tr_seri)
891
892	elseif (tr_scheme.eq.2) then
893	! Case 2: surface emission
894	! For the moment, inspired from Mars version
895	! However, the variable 'source' could be used in physiq
896	! so the call to phytrac_emiss could be to initialise it.
897
898	call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
899	I debut,lafin,nqmax,
900	I nlon,nlev,dtime,paprs,
901	I latitude_deg,longitude_deg,
902	O tr_seri)
903
904	elseif (tr_scheme.eq.3) then ! identical to ok_chem.or.ok_cloud
905	! Case 3: Full chemistry and/or clouds
906
907	call phytrac_chimie(
908	I debut,
909	I gmtime,
910	I nqmax,
911	I klon,
912	I latitude_deg,
913	I longitude_deg,
914	I nlev,
915	I dtime,
916	I t_seri,
917	I pplay,
918	O tr_seri)
919
920	c CALL WriteField_phy('Pression',pplay,nlev)
921	c CALL WriteField_phy('PressionBnd',paprs,nlev+1)
922	c CALL WriteField_phy('Temp',t_seri,nlev)
923	c IF (ok_cloud) THEN
924	c CALL WriteField_phy('NBRTOT',NBRTOT,nlev)
925	c ENDIF
926	c CALL WriteField_phy('SAl',tr_seri(:,:,i_h2so4liq),nlev)
927	c CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev)
928
929	C === SEDIMENTATION ===
930
931	if (ok_sedim) then
932
933	c !! Depend on cl_scheme !!
934
935	if (cl_scheme.eq.1) then
936	c ================
937
938	CALL new_cloud_sedim(
939	I klon,
940	I nlev,
941	I dtime,
942	I pplay,
943	I paprs,
944	I t_seri,
945	I tr_seri,
946	O d_tr_sed(:,:,1:2),
947	O d_tr_ssed,
948	I nqmax,
949	O Fsedim)
950
951	DO k = 1, klev
952	DO i = 1, klon
953
954	c--------------------
955	c Ce test est necessaire pour eviter Xliq=NaN
956	! IF (ieee_is_nan(d_tr_sed(i,k,1)).OR.
957	! & ieee_is_nan(d_tr_sed(i,k,2))) THEN
958	IF ((d_tr_sed(i,k,1).ne.d_tr_sed(i,k,1)).OR.
959	& (d_tr_sed(i,k,2).ne.d_tr_sed(i,k,2))) THEN
960	PRINT*,'sedim NaN PROBLEM'
961	PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k)
962	PRINT*,'lat-lon',i,'level',k,'dtime',dtime
963	PRINT*,'F_sed',Fsedim(i,k)
964	PRINT*,'==============================================='
965	d_tr_sed(i,k,:)=0.
966	ENDIF
967	c--------------------
968
969	tr_seri(i,k,i_h2so4liq) = tr_seri(i,k,i_h2so4liq)+
970	& d_tr_sed(i,k,1)
971	tr_seri(i,k,i_h2oliq) = tr_seri(i,k,i_h2oliq)+
972	& d_tr_sed(i,k,2)
973	d_tr_sed(i,k,:) = d_tr_sed(i,k,:) / dtime
974	Fsedim(i,k) = Fsedim(i,k) / dtime
975
976	ENDDO
977	ENDDO
978
979	Fsedim(:,klev+1) = 0.
980
981	elseif (cl_scheme.eq.2) then
982	c ====================
983
984	d_tr_sed(:,:,:) = 0.D0
985
986	DO i=1, klon
987
988	c Mode 1
989	m0_mode1(:,1)=tr_seri(i,:,i_m0_mode1drop)
990	m0_mode1(:,2)=tr_seri(i,:,i_m0_mode1ccn)
991	m3_mode1(:,1)=tr_seri(i,:,i_m3_mode1sa)
992	m3_mode1(:,2)=tr_seri(i,:,i_m3_mode1w)
993	m3_mode1(:,3)=tr_seri(i,:,i_m3_mode1ccn)
994
995	call drop_sedimentation(dtime,klev,m0_mode1,m3_mode1,
996	& paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:),
997	& 1,
998	& d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed)
999
1000	d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop)
1001	& + d_drop_sed
1002	d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn)
1003	& + d_ccn_sed(:,1)
1004	d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn)
1005	& + d_ccn_sed(:,2)
1006	d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa)
1007	& + d_liq_sed(:,1)
1008	d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w)
1009	& + d_liq_sed(:,2)
1010
1011	c Mode 2
1012	m0_mode2(:,1)=tr_seri(i,:,i_m0_mode2drop)
1013	m0_mode2(:,2)=tr_seri(i,:,i_m0_mode2ccn)
1014	m3_mode2(:,1)=tr_seri(i,:,i_m3_mode2sa)
1015	m3_mode2(:,2)=tr_seri(i,:,i_m3_mode2w)
1016	m3_mode2(:,3)=tr_seri(i,:,i_m3_mode2ccn)
1017
1018	call drop_sedimentation(dtime,klev,m0_mode2,m3_mode2,
1019	& paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:),
1020	& 2,
1021	& d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed)
1022
1023	d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop)
1024	& + d_drop_sed
1025	d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn)
1026	& + d_ccn_sed(:,1)
1027	d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn)
1028	& + d_ccn_sed(:,2)
1029	d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa)
1030	& + d_liq_sed(:,1)
1031	d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w)
1032	& + d_liq_sed(:,2)
1033
1034	c Aer
1035	call aer_sedimentation(dtime,klev,tr_seri(i,:,i_m0_aer),
1036	& tr_seri(i,:,i_m3_aer),paprs(i,:),
1037	& zzlev(i,:),zzlay(i,:),t_seri(i,:),
1038	& d_tr_sed(i,:,i_m0_aer),d_tr_sed(i,:,i_m3_aer),aer_flux)
1039
1040	END DO
1041
1042	DO iq=nqmax-nmicro+1,nqmax
1043	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_sed(:,:,iq)
1044	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq) / dtime
1045	END DO
1046
1047	endif
1048	c ====================
1049
1050	C === FIN SEDIMENTATION ===
1051
1052	endif ! ok_sedim
1053
1054	endif ! tr_scheme
1055	endif ! iflag_trac
1056
1057	c
1058	c====================================================================
1059	c Appeler la diffusion verticale (programme de couche limite)
1060	c====================================================================
1061
1062	c-------------------------------
1063	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
1064	c l'equilibre radiatif du sol
1065	if (1.eq.0) then
1066	if (debut) then
1067	print*,"ATTENTION, CLMAIN SHUNTEE..."
1068	endif
1069
1070	DO i = 1, klon
1071	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
1072	fder(i) = 0.0e0
1073	dlw(i) = 0.0e0
1074	ENDDO
1075
1076	c Incrementer la temperature du sol
1077	c
1078	DO i = 1, klon
1079	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
1080	ftsol(i) = ftsol(i) + d_ts(i)
1081	do j=1,nsoilmx
1082	ftsoil(i,j)=ftsol(i)
1083	enddo
1084	ENDDO
1085
1086	c-------------------------------
1087	else
1088	c-------------------------------
1089
1090	fder = dlw
1091
1092	! ADAPTATION GCM POUR CP(T)
1093
1094	CALL clmain(dtime,itap,
1095	e t_seri,u_seri,v_seri,
1096	e rmu0,
1097	e ftsol,
1098	$ ftsoil,
1099	$ paprs,pplay,ppk,radsol,falbe,
1100	e solsw, sollw, sollwdown, fder,
1101	e longitude_deg, latitude_deg, dx, dy,
1102	e debut, lafin,
1103	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1104	s fluxt,fluxu,fluxv,cdragh,cdragm,
1105	s dsens,
1106	s ycoefh,yu1,yv1)
1107
1108	CXXX Incrementation des flux
1109	DO i = 1, klon
1110	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
1111	fder(i) = dlw(i) + dsens(i)
1112	ENDDO
1113	CXXX
1114
1115	DO k = 1, klev
1116	DO i = 1, klon
1117	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1118	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
1119	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1120	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
1121	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1122	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
1123	ENDDO
1124	ENDDO
1125
1126	C TRACEURS
1127
1128	if (iflag_trac.eq.1) then
1129	DO k = 1, klev
1130	DO i = 1, klon
1131	delp(i,k) = paprs(i,k)-paprs(i,k+1)
1132	ENDDO
1133	ENDDO
1134
1135	DO iq=1, nqmax
1136
1137	CALL cltrac(dtime,ycoefh,t_seri,
1138	s tr_seri(:,:,iq),source(:,iq),
1139	e paprs, pplay,delp,
1140	s d_tr_vdf(:,:,iq))
1141
1142	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
1143	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
1144
1145	ENDDO !nqmax
1146
1147	endif
1148
1149	IF (if_ebil.ge.2) THEN
1150	ztit='after clmain'
1151	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
1152	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1153	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1154	call diagphy(cell_area,ztit,ip_ebil
1155	e , zero_v, zero_v, zero_v, zero_v, sens
1156	e , zero_v, zero_v, zero_v, ztsol
1157	e , d_h_vcol, d_qt, d_ec
1158	s , fs_bound, fq_bound )
1159	END IF
1160	C
1161	c
1162	c Incrementer la temperature du sol
1163	c
1164	DO i = 1, klon
1165	ftsol(i) = ftsol(i) + d_ts(i)
1166	ENDDO
1167
1168	c Calculer la derive du flux infrarouge
1169	c
1170	DO i = 1, klon
1171	dlw(i) = - 4.0RSIGMAftsol(i)**3
1172	ENDDO
1173
1174	c-------------------------------
1175	endif ! fin du VENUS TEST
1176
1177	! tests: output tendencies
1178	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
1179	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
1180	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
1181	! call writefield_phy('physiq_d_ts',d_ts,1)
1182
1183	c
1184	c Appeler l'ajustement sec
1185	c
1186	c===================================================================
1187	c Convection seche
1188	c===================================================================
1189	c
1190	d_t_ajs(:,:)=0.
1191	d_u_ajs(:,:)=0.
1192	d_v_ajs(:,:)=0.
1193	d_tr_ajs(:,:,:)=0.
1194	c
1195	IF(prt_level>9)WRITE(lunout,*)
1196	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
1197	s ,iflag_ajs
1198
1199	if(iflag_ajs.eq.0) then
1200	c Rien
1201	c ====
1202	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1203
1204	else if(iflag_ajs.eq.1) then
1205
1206	c Ajustement sec
1207	c ==============
1208	IF(prt_level>9)WRITE(lunout,*)'ajsec'
1209
1210	! ADAPTATION GCM POUR CP(T)
1211	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
1212	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
1213
1214	! ADAPTATION GCM POUR CP(T)
1215	do i=1,klon
1216	flux_ajs(i,1) = 0.0
1217	do j=2,klev
1218	flux_ajs(i,j) = flux_ajs(i,j-1)
1219	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
1220	. *(paprs(i,j-1)-paprs(i,j))
1221	enddo
1222	enddo
1223
1224	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1225	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1226	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1227	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1228	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1229	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1230
1231	if (iflag_trac.eq.1) then
1232	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1233	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1234	endif
1235	endif
1236
1237	! tests: output tendencies
1238	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
1239	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
1240	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
1241	c
1242	IF (if_ebil.ge.2) THEN
1243	ztit='after dry_adjust'
1244	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1245	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1246	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1247	call diagphy(cell_area,ztit,ip_ebil
1248	e , zero_v, zero_v, zero_v, zero_v, sens
1249	e , zero_v, zero_v, zero_v, ztsol
1250	e , d_h_vcol, d_qt, d_ec
1251	s , fs_bound, fq_bound )
1252	END IF
1253
1254
1255	c====================================================================
1256	c RAYONNEMENT
1257	c====================================================================
1258
1259	c------------------------------------
1260	c . Compute radiative tendencies :
1261	c------------------------------------
1262	c====================================================================
1263	IF (MOD(itaprad,radpas).EQ.0) THEN
1264	c====================================================================
1265
1266	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1267	c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
1268
1269
1270	c------------------------------------
1271	c . Compute mean mass, cp and R :
1272	c------------------------------------
1273
1274	if(callthermos) then
1275	call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax)
1276
1277	endif
1278
1279
1280	cc!!! ADD key callhedin
1281
1282	IF(callnlte.or.callthermos) THEN
1283	call compo_hedin83_mod(pplay,rmu0,
1284	& co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1285
1286	IF(ok_chem) then
1287
1288	CC !! GG : Using only mayor species tracers abundances to compute NLTE heating/cooling
1289
1290	CC Conversion [mmr] ---> [vmr]
1291
1292	co2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co2)*
1293	& mmean(1:nlon,1:nlev)/M_tr(i_co2)
1294	covmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co)*
1295	& mmean(1:nlon,1:nlev)/M_tr(i_co)
1296	ovmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_o)*
1297	& mmean(1:nlon,1:nlev)/M_tr(i_o)
1298	n2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_n2)*
1299	& mmean(1:nlon,1:nlev)/M_tr(i_n2)
1300
1301	ENDIF
1302
1303	ENDIF
1304
1305	c
1306	c NLTE cooling from CO2 emission
1307	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1308
1309	IF(callnlte) THEN
1310	if(nltemodel.eq.0.or.nltemodel.eq.1) then
1311	CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
1312	$ tr_seri, d_t_nlte)
1313	else if(nltemodel.eq.2) then
1314	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
1315	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1316	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
1317	if(ierr_nlte.gt.0) then
1318	write(,)
1319	$ 'WARNING: nlte_tcool output with error message',
1320	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
1321	write(,)'I will continue anyway'
1322	endif
1323
1324	endif
1325
1326	ELSE
1327
1328	d_t_nlte(:,:)=0.
1329
1330	ENDIF
1331
1332	c Find number of layers for LTE radiation calculations
1333
1334	IF(callnlte .or. callnirco2)
1335	$ CALL nlthermeq(klon, klev, paprs, pplay)
1336
1337	c
1338	c LTE radiative transfert / solar / IR matrix
1339	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1340	CALL radlwsw
1341	e (dist, rmu0, fract, zzlev,
1342	e paprs, pplay,ftsol, t_seri)
1343
1344
1345	c CO2 near infrared absorption
1346	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1347
1348	d_t_nirco2(:,:)=0.
1349	if (callnirco2) then
1350	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
1351	. rmu0, fract, d_t_nirco2)
1352	endif
1353
1354
1355	c Net atmospheric radiative heating rate (K.s-1)
1356	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1357
1358	IF(callnlte.or.callnirco2) THEN
1359	CALL blendrad(klon, klev, pplay,heat,
1360	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
1361	ELSE
1362	dtsw(:,:)=heat(:,:)
1363	dtlw(:,:)=-1*cool(:,:)
1364	ENDIF
1365
1366	DO k=1,klev
1367	DO i=1,klon
1368	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
1369	ENDDO
1370	ENDDO
1371
1372
1373	cc---------------------------------------------
1374
1375	c EUV heating rate (K.s-1)
1376	c ~~~~~~~~~~~~~~~~~~~~~~~~
1377
1378	d_t_euv(:,:)=0.
1379
1380	IF (callthermos) THEN
1381
1382	c call euvheat(klon, klev,t_seri,paprs,pplay,zzlay,
1383	c $ rmu0,pdtphys,gmtime,rjourvrai, co2vmr_gcm, n2vmr_gcm,
1384	c $ covmr_gcm, ovmr_gcm,d_t_euv )
1385	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
1386	$ rmu0,pdtphys,gmtime,rjourvrai,
1387	$ tr_seri, d_tr, d_t_euv )
1388
1389	DO k=1,klev
1390	DO ig=1,klon
1391	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
1392
1393	ENDDO
1394	ENDDO
1395
1396	ENDIF ! callthermos
1397
1398	c====================================================================
1399	itaprad = 0
1400	ENDIF ! radpas
1401	c====================================================================
1402	c
1403	c Ajouter la tendance des rayonnements (tous les pas)
1404	c
1405	DO k = 1, klev
1406	DO i = 1, klon
1407	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
1408	ENDDO
1409	ENDDO
1410
1411	! CONDUCTION and MOLECULAR VISCOSITY
1412	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1413
1414	d_t_conduc(:,:)=0.
1415	d_u_molvis(:,:)=0.
1416	d_v_molvis(:,:)=0.
1417
1418	IF (callthermos) THEN
1419
1420	tsurf(:)=t_seri(:,1)
1421	call conduction(klon, klev,pdtphys,
1422	$ pplay,paprs,t_seri,
1423	$ tsurf,zzlev,zzlay,d_t_conduc)
1424
1425	call molvis(klon, klev,pdtphys,
1426	$ pplay,paprs,t_seri,
1427	$ u,tsurf,zzlev,zzlay,d_u_molvis)
1428
1429	call molvis(klon, klev, pdtphys,
1430	$ pplay,paprs,t_seri,
1431	$ v,tsurf,zzlev,zzlay,d_v_molvis)
1432
1433	DO k=1,klev
1434	DO ig=1,klon
1435	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
1436	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
1437	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
1438	ENDDO
1439	ENDDO
1440	ENDIF
1441
1442
1443	! -- MOLECULAR DIFFUSION ---
1444
1445	d_q_moldif(:,:,:)=0
1446
1447	IF (callthermos .and. ok_chem) THEN
1448
1449	call moldiff_red(klon, klev, nqmax,
1450	& pplay,paprs,t_seri, tr_seri, pdtphys,
1451	& zzlay,d_t_euv,d_t_conduc,d_q_moldif)
1452
1453
1454	! --- update tendencies tracers ---
1455
1456	DO iq = 1, nqmax
1457	DO k=1,klev
1458	DO ig=1,klon
1459	tr_seri(ig,k,iq)= tr_seri(ig,k,iq)+
1460	& d_q_moldif(ig,k,iq)*dtime ! [Kg/kg]?
1461	ENDDO
1462	ENDDO
1463	ENDDO
1464
1465
1466	ENDIF ! callthermos & ok_chem
1467
1468	c====================================================================
1469	! tests: output tendencies
1470	! call writefield_phy('physiq_dtrad',dtrad,klev)
1471
1472	IF (if_ebil.ge.2) THEN
1473	ztit='after rad'
1474	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1475	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1476	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1477	call diagphy(cell_area,ztit,ip_ebil
1478	e , topsw, toplw, solsw, sollw, zero_v
1479	e , zero_v, zero_v, zero_v, ztsol
1480	e , d_h_vcol, d_qt, d_ec
1481	s , fs_bound, fq_bound )
1482	END IF
1483	c
1484
1485	c====================================================================
1486	c Calcul des gravity waves FLOTT
1487	c====================================================================
1488	c
1489	if (ok_orodr.or.ok_gw_nonoro) then
1490
1491	c CALCUL DE N2
1492	c UTILISE LA RELATION ENTRE N2 ET STABILITE
1493	c N2 = RG/T (dT/dz+RG/cp(T))
1494	c ET DONC EN N'UTILISE QUE LA TEMPERATURE, PAS teta.
1495
1496	do i=1,klon
1497	do k=2,klev
1498	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1499	enddo
1500	enddo
1501	do i=1,klon
1502	do k=2,klev
1503	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1504	zdtlev(i,k) = t_seri(i,k)-t_seri(i,k-1)
1505	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
1506	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
1507	. + RG/cpdet(ztlev(i,k)) )
1508	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1509	enddo
1510	zn2(i,1) = 1.e-12 ! securite
1511	enddo
1512
1513	endif
1514
1515	c ----------------------------ORODRAG
1516	IF (ok_orodr) THEN
1517	c
1518	c selection des points pour lesquels le shema est actif:
1519	igwd=0
1520	DO i=1,klon
1521	itest(i)=0
1522	c IF ((zstd(i).gt.10.0)) THEN
1523	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1524	itest(i)=1
1525	igwd=igwd+1
1526	idx(igwd)=i
1527	ENDIF
1528	ENDDO
1529	c igwdim=MAX(1,igwd)
1530	c
1531	c A ADAPTER POUR VENUS!!!
1532	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
1533	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1534	e igwd,idx,itest,
1535	e t_seri, u_seri, v_seri,
1536	s zulow, zvlow, zustrdr, zvstrdr,
1537	s d_t_oro, d_u_oro, d_v_oro)
1538
1539	c print*,"d_u_oro=",d_u_oro(klon/2,:)
1540	c ajout des tendances
1541	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1542	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1543	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1544	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1545	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1546	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1547	c
1548	ELSE
1549	d_t_oro = 0.
1550	d_u_oro = 0.
1551	d_v_oro = 0.
1552	zustrdr = 0.
1553	zvstrdr = 0.
1554	c
1555	ENDIF ! fin de test sur ok_orodr
1556	c
1557	! tests: output tendencies
1558	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
1559	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
1560	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
1561
1562	c ----------------------------OROLIFT
1563	IF (ok_orolf) THEN
1564	print*,"ok_orolf NOT IMPLEMENTED !"
1565	stop
1566	c
1567	c selection des points pour lesquels le shema est actif:
1568	igwd=0
1569	DO i=1,klon
1570	itest(i)=0
1571	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1572	itest(i)=1
1573	igwd=igwd+1
1574	idx(igwd)=i
1575	ENDIF
1576	ENDDO
1577	c igwdim=MAX(1,igwd)
1578	c
1579	c A ADAPTER POUR VENUS!!!
1580	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
1581	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1582	c e igwd,idx,itest,
1583	c e t_seri, u_seri, v_seri,
1584	c s zulow, zvlow, zustrli, zvstrli,
1585	c s d_t_lif, d_u_lif, d_v_lif )
1586
1587	c
1588	c ajout des tendances
1589	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1590	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1591	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1592	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1593	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1594	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1595	c
1596	ELSE
1597	d_t_lif = 0.
1598	d_u_lif = 0.
1599	d_v_lif = 0.
1600	zustrli = 0.
1601	zvstrli = 0.
1602	c
1603	ENDIF ! fin de test sur ok_orolf
1604
1605	c ---------------------------- NON-ORO GRAVITY WAVES
1606	IF(ok_gw_nonoro) then
1607
1608	call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1609	e t_seri, u_seri, v_seri, plevmoy,
1610	o zustrhi,zvstrhi,
1611	o d_t_hin, d_u_hin, d_v_hin)
1612
1613	c ajout des tendances
1614
1615	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1616	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1617	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1618	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1619	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1620	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1621
1622	ELSE
1623	d_t_hin = 0.
1624	d_u_hin = 0.
1625	d_v_hin = 0.
1626	zustrhi = 0.
1627	zvstrhi = 0.
1628
1629	ENDIF ! fin de test sur ok_gw_nonoro
1630
1631	! tests: output tendencies
1632	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
1633	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
1634	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
1635
1636	c====================================================================
1637	c Transport de ballons
1638	c====================================================================
1639	if (ballons.eq.1) then
1640	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
1641	& latitude_deg,longitude_deg,
1642	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1643	C t,pplay,u,v,zphi) ! alt above planet average radius
1644	endif !ballons
1645
1646	c====================================================================
1647	c Bilan de mmt angulaire
1648	c====================================================================
1649	if (bilansmc.eq.1) then
1650	CMODDEB FLOTT
1651	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1652	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1653
1654	DO i = 1, klon
1655	zustrph(i)=0.
1656	zvstrph(i)=0.
1657	zustrcl(i)=0.
1658	zvstrcl(i)=0.
1659	ENDDO
1660	DO k = 1, klev
1661	DO i = 1, klon
1662	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1663	c (paprs(i,k)-paprs(i,k+1))/rg
1664	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1665	c (paprs(i,k)-paprs(i,k+1))/rg
1666	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1667	c (paprs(i,k)-paprs(i,k+1))/rg
1668	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1669	c (paprs(i,k)-paprs(i,k+1))/rg
1670	ENDDO
1671	ENDDO
1672
1673	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
1674	C ra,rg,romega,
1675	C latitude_deg,longitude_deg,pphis,
1676	C zustrdr,zustrli,zustrcl,
1677	C zvstrdr,zvstrli,zvstrcl,
1678	C paprs,u,v)
1679
1680	CCMODFIN FLOTT
1681	endif !bilansmc
1682
1683	c====================================================================
1684	c====================================================================
1685	c Calculer le transport de l'eau et de l'energie (diagnostique)
1686	c
1687	c A REVOIR POUR VENUS...
1688	c
1689	c CALL transp (paprs,ftsol,
1690	c e t_seri, q_seri, u_seri, v_seri, zphi,
1691	c s ve, vq, ue, uq)
1692	c
1693	c====================================================================
1694	c+jld ec_conser
1695	DO k = 1, klev
1696	DO i = 1, klon
1697	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1698	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1699	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1700	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1701	END DO
1702	END DO
1703	do i=1,klon
1704	flux_ec(i,1) = 0.0
1705	do j=2,klev
1706	flux_ec(i,j) = flux_ec(i,j-1)
1707	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
1708	enddo
1709	enddo
1710
1711	c-jld ec_conser
1712	c====================================================================
1713	IF (if_ebil.ge.1) THEN
1714	ztit='after physic'
1715	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
1716	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1717	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1718	C Comme les tendances de la physique sont ajoute dans la dynamique,
1719	C on devrait avoir que la variation d'entalpie par la dynamique
1720	C est egale a la variation de la physique au pas de temps precedent.
1721	C Donc la somme de ces 2 variations devrait etre nulle.
1722	call diagphy(cell_area,ztit,ip_ebil
1723	e , topsw, toplw, solsw, sollw, sens
1724	e , zero_v, zero_v, zero_v, ztsol
1725	e , d_h_vcol, d_qt, d_ec
1726	s , fs_bound, fq_bound )
1727	C
1728	d_h_vcol_phy=d_h_vcol
1729	C
1730	END IF
1731	C
1732	c=======================================================================
1733	c SORTIES
1734	c=======================================================================
1735
1736	c Convertir les incrementations en tendances
1737	c
1738	DO k = 1, klev
1739	DO i = 1, klon
1740	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1741	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1742	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1743	ENDDO
1744	ENDDO
1745	c
1746	DO iq = 1, nqmax
1747	DO k = 1, klev
1748	DO i = 1, klon
1749	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1750	ENDDO
1751	ENDDO
1752	ENDDO
1753
1754	c------------------------
1755	c Calcul moment cinetique
1756	c------------------------
1757	c TEST VENUS...
1758	c mangtot = 0.0
1759	c DO k = 1, klev
1760	c DO i = 1, klon
1761	c mang(i,k) = RA*cos(latitude(i))
1762	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
1763	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
1764	c mangtot=mangtot+mang(i,k)
1765	c ENDDO
1766	c ENDDO
1767	c print*,"Moment cinetique total = ",mangtot
1768	c
1769	c------------------------
1770	c
1771	c Sauvegarder les valeurs de t et u a la fin de la physique:
1772	c
1773	DO k = 1, klev
1774	DO i = 1, klon
1775	u_ancien(i,k) = u_seri(i,k)
1776	t_ancien(i,k) = t_seri(i,k)
1777	ENDDO
1778	ENDDO
1779	c
1780	c=============================================================
1781	c Ecriture des sorties
1782	c=============================================================
1783
1784	#ifdef CPP_IOIPSL
1785
1786	#ifdef histhf
1787	#include "write_histhf.h"
1788	#endif
1789
1790	#ifdef histday
1791	#include "write_histday.h"
1792	#endif
1793
1794	#ifdef histmth
1795	#include "write_histmth.h"
1796	#endif
1797
1798	#ifdef histins
1799	#include "write_histins.h"
1800	#endif
1801
1802	#endif
1803
1804	! XIOS outputs
1805	! This can be done ANYWHERE in the physics routines !
1806
1807	#ifdef CPP_XIOS
1808	! Send fields to XIOS: (NB these fields must also be defined as
1809	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
1810
1811	! 2D fields
1812
1813	CALL send_xios_field("phis",pphis)
1814	cell_area_out(:)=cell_area(:)
1815	if (is_north_pole_phy) cell_area_out(1)=cell_area(1)/nbp_lon
1816	if (is_south_pole_phy) cell_area_out(klon)=cell_area(klon)/nbp_lon
1817	CALL send_xios_field("aire",cell_area_out)
1818	CALL send_xios_field("tsol",ftsol)
1819	CALL send_xios_field("psol",paprs(:,1))
1820	CALL send_xios_field("cdragh",cdragh)
1821	CALL send_xios_field("cdragm",cdragm)
1822
1823	CALL send_xios_field("tops",topsw)
1824	CALL send_xios_field("topl",toplw)
1825	CALL send_xios_field("sols",solsw)
1826	CALL send_xios_field("soll",sollw)
1827
1828	! 3D fields
1829
1830	CALL send_xios_field("temp",t_seri)
1831	CALL send_xios_field("pres",pplay)
1832	CALL send_xios_field("geop",zphi)
1833	CALL send_xios_field("vitu",u_seri)
1834	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1835	CALL send_xios_field("vitv",-1.*v_seri)
1836	CALL send_xios_field("vitw",omega)
1837	CALL send_xios_field("Kz",ycoefh)
1838	CALL send_xios_field("mmean",mmean)
1839	CALL send_xios_field("rho",rho)
1840
1841	CALL send_xios_field("dudyn",d_u_dyn)
1842	CALL send_xios_field("duvdf",d_u_vdf)
1843	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1844	CALL send_xios_field("dvvdf",-1.*d_v_vdf)
1845	CALL send_xios_field("duajs",d_u_ajs)
1846	CALL send_xios_field("dugwo",d_u_oro)
1847	CALL send_xios_field("dugwno",d_u_hin)
1848	CALL send_xios_field("dumolvis",d_u_molvis)
1849	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1850	CALL send_xios_field("dvmolvis",-1.*d_v_molvis)
1851	CALL send_xios_field("dtdyn",d_t_dyn)
1852	CALL send_xios_field("dtphy",d_t)
1853	CALL send_xios_field("dtvdf",d_t_vdf)
1854	CALL send_xios_field("dtajs",d_t_ajs)
1855	CALL send_xios_field("dtswr",dtsw)
1856	CALL send_xios_field("dtswrNLTE",d_t_nirco2)
1857	CALL send_xios_field("dtswrLTE",heat)
1858	CALL send_xios_field("dtlwr",dtlw)
1859	CALL send_xios_field("dtlwrNLTE",d_t_nlte)
1860	CALL send_xios_field("dtlwrLTE",-1.*cool)
1861	CALL send_xios_field("dteuv",d_t_euv)
1862	CALL send_xios_field("dtcond",d_t_conduc)
1863	CALL send_xios_field("dtec",d_t_ec)
1864
1865	CALL send_xios_field("SWnet",swnet(:,1:klev))
1866	CALL send_xios_field("LWnet",lwnet(:,1:klev))
1867	CALL send_xios_field("fluxvdf",fluxt)
1868	CALL send_xios_field("fluxdyn",flux_dyn)
1869	CALL send_xios_field("fluxajs",flux_ajs)
1870	CALL send_xios_field("fluxec",flux_ec)
1871
1872	c When using tracers
1873	IF (iflag_trac.eq.1) THEN
1874	c photochemical compounds !!!outputs in [vmr]
1875	DO iq=1,nqmax-nmicro
1876	CALL send_xios_field(tname(iq),qx(:,:,iq)*mmean(:,:)/M_tr(iq))
1877	ENDDO
1878
1879	IF ((tr_scheme.eq.3).and.(ok_sedim).and.(cl_scheme.eq.1)) THEN
1880	CALL send_xios_field("Fsedim",Fsedim(:,1:klev))
1881	ENDIF
1882	ENDIF
1883
1884	IF (callthermos .and. ok_chem) THEN
1885	CALL send_xios_field("d_qmoldifCO2",d_q_moldif(:,:,i_co2))
1886	CALL send_xios_field("d_qmoldifO3p",d_q_moldif(:,:,i_o))
1887	CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
1888	ENDIF
1889
1890	if (lafin.and.is_omp_master) then
1891	write(,) "physiq: call xios_context_finalize"
1892	call xios_context_finalize
1893	endif
1894
1895	#endif
1896
1897	c====================================================================
1898	c Si c'est la fin, il faut conserver l'etat de redemarrage
1899	c====================================================================
1900	c
1901	IF (lafin) THEN
1902	itau_phy = itau_phy + itap
1903	CALL phyredem ("restartphy.nc")
1904
1905	c--------------FLOTT
1906	CMODEB LOTT
1907	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
1908	C DU BILAN DE MOMENT ANGULAIRE.
1909	if (bilansmc.eq.1) then
1910	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
1911	close(27)
1912	close(28)
1913	endif !bilansmc
1914	CMODFIN
1915	c-------------
1916	c--------------SLEBONNOIS
1917	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
1918	C DES BALLONS
1919	if (ballons.eq.1) then
1920	write(,)'Fermeture des ballons*.out'
1921	close(30)
1922	close(31)
1923	close(32)
1924	close(33)
1925	close(34)
1926	endif !ballons
1927	c-------------
1928	ENDIF
1929
1930	END SUBROUTINE physiq
1931
1932	END MODULE physiq_mod
1933

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

Download in other formats:

Original Format