Context Navigation

physiq_mod.F @ 1676

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

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

Download in other formats:

Original Format