Context Navigation

physiq_mod.F @ 2193

Last change on this file since 2193 was 2193, checked in by flefevre, 5 years ago
mise en place du supercycling du rayonnement et de la chimie
File size: 64.7 KB

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

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

Download in other formats:

Original Format