Context Navigation

physiq_mod.F @ 2194

Last change on this file since 2194 was 2194, checked in by slebonnois, 5 years ago
SL: ajustement commit precedent pour la microphysique de Sabrina
File size: 65.1 KB

Line
1	!
2	! $Id: $
3	!
4	MODULE physiq_mod
5
6	IMPLICIT NONE
7
8	CONTAINS
9
10	SUBROUTINE physiq (nlon,nlev,nqmax,
11	. debut,lafin,rjourvrai,gmtime,pdtphys,
12	. paprs,pplay,ppk,pphi,pphis,presnivs,
13	. u,v,t,qx,
14	. flxmw,
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	d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime
1117	end do
1118	#endif
1119	end if ! cl_scheme
1120	end if ! ok_sedim
1121	end if ! mod(itap,chempas) <------- end of chemistry supercycling
1122
1123	! update tracers (chemistry)
1124
1125	do iq = 1, nqmax - nmicro
1126	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_chem(:,:,iq)*dtime
1127	end do
1128
1129	! update tracers (sedimentation)
1130
1131	if (ok_sedim) then
1132	if (cl_scheme == 1) then
1133	tr_seri(:,:,i_h2so4liq) = tr_seri(:,:,i_h2so4liq)
1134	$ + d_tr_sed(:,:,1)*dtime
1135	tr_seri(:,:,i_h2oliq) = tr_seri(:,:,i_h2oliq)
1136	$ + d_tr_sed(:,:,2)*dtime
1137	else if (cl_scheme == 2) then
1138	do iq = nqmax-nmicro+1,nqmax
1139	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_sed(:,:,iq)*dtime
1140	end do
1141	end if ! cl_scheme
1142	end if ! ok_sedim
1143	!====================================================================
1144	! End Case 3: Full chemistry and/or clouds.
1145	!====================================================================
1146
1147	end if ! tr_scheme
1148	end if ! iflag_trac
1149
1150	c====================================================================
1151	c Appeler la diffusion verticale (programme de couche limite)
1152	c====================================================================
1153
1154	c-------------------------------
1155	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
1156	c l'equilibre radiatif du sol
1157	if (.not. ok_clmain) then
1158	if (debut) then
1159	print*,"ATTENTION, CLMAIN SHUNTEE..."
1160	endif
1161
1162	DO i = 1, klon
1163	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
1164	fder(i) = 0.0e0
1165	dlw(i) = 0.0e0
1166	ENDDO
1167
1168	c Incrementer la temperature du sol
1169	c
1170	DO i = 1, klon
1171	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
1172	ftsol(i) = ftsol(i) + d_ts(i)
1173	do j=1,nsoilmx
1174	ftsoil(i,j)=ftsol(i)
1175	enddo
1176	ENDDO
1177
1178	c-------------------------------
1179	else
1180	c-------------------------------
1181
1182	fder = dlw
1183
1184	! ADAPTATION GCM POUR CP(T)
1185
1186	if (physideal) then
1187	CALL clmain_ideal(dtime,itap,
1188	e t_seri,u_seri,v_seri,
1189	e rmu0,
1190	e ftsol,
1191	$ ftsoil,
1192	$ paprs,pplay,ppk,radsol,falbe,
1193	e solsw, sollw, sollwdown, fder,
1194	e longitude_deg, latitude_deg, dx, dy,
1195	e debut, lafin,
1196	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1197	s fluxt,fluxu,fluxv,cdragh,cdragm,
1198	s dsens,
1199	s ycoefh,yu1,yv1)
1200	else
1201	CALL clmain(dtime,itap,
1202	e t_seri,u_seri,v_seri,
1203	e rmu0,
1204	e ftsol,
1205	$ ftsoil,
1206	$ paprs,pplay,ppk,radsol,falbe,
1207	e solsw, sollw, sollwdown, fder,
1208	e longitude_deg, latitude_deg, dx, dy,
1209	e debut, lafin,
1210	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
1211	s fluxt,fluxu,fluxv,cdragh,cdragm,
1212	s dsens,
1213	s ycoefh,yu1,yv1)
1214	endif
1215
1216	CXXX Incrementation des flux
1217	DO i = 1, klon
1218	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
1219	fder(i) = dlw(i) + dsens(i)
1220	ENDDO
1221	CXXX
1222	IF (.not. turb_resolved) then !True only for LES
1223	DO k = 1, klev
1224	DO i = 1, klon
1225	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1226	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
1227	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1228	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
1229	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1230	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
1231	ENDDO
1232	ENDDO
1233	ENDIF
1234	C TRACEURS
1235
1236	if (iflag_trac.eq.1) then
1237	DO k = 1, klev
1238	DO i = 1, klon
1239	delp(i,k) = paprs(i,k)-paprs(i,k+1)
1240	ENDDO
1241	ENDDO
1242
1243	DO iq=1, nqmax
1244
1245	CALL cltrac(dtime,ycoefh,t_seri,
1246	s tr_seri(:,:,iq),source(:,iq),
1247	e paprs, pplay,delp,
1248	s d_tr_vdf(:,:,iq))
1249
1250	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
1251	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
1252
1253	ENDDO !nqmax
1254
1255	endif
1256
1257	IF (if_ebil.ge.2) THEN
1258	ztit='after clmain'
1259	CALL diagetpq(cell_area,ztit,ip_ebil,2,1,dtime
1260	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1261	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1262	call diagphy(cell_area,ztit,ip_ebil
1263	e , zero_v, zero_v, zero_v, zero_v, sens
1264	e , zero_v, zero_v, zero_v, ztsol
1265	e , d_h_vcol, d_qt, d_ec
1266	s , fs_bound, fq_bound )
1267	END IF
1268	C
1269	c
1270	c Incrementer la temperature du sol
1271	c
1272	DO i = 1, klon
1273	ftsol(i) = ftsol(i) + d_ts(i)
1274	ENDDO
1275
1276	c Calculer la derive du flux infrarouge
1277	c
1278	DO i = 1, klon
1279	dlw(i) = - 4.0RSIGMAftsol(i)**3
1280	ENDDO
1281
1282	c-------------------------------
1283	endif ! fin du VENUS TEST
1284
1285	! tests: output tendencies
1286	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
1287	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
1288	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
1289	! call writefield_phy('physiq_d_ts',d_ts,1)
1290
1291	c
1292	c Appeler l'ajustement sec
1293	c
1294	c===================================================================
1295	c Convection seche
1296	c===================================================================
1297	c
1298	d_t_ajs(:,:)=0.
1299	d_u_ajs(:,:)=0.
1300	d_v_ajs(:,:)=0.
1301	d_tr_ajs(:,:,:)=0.
1302	c
1303	IF(prt_level>9)WRITE(lunout,*)
1304	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
1305	s ,iflag_ajs
1306
1307	if(iflag_ajs.eq.0) then
1308	c Rien
1309	c ====
1310	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1311
1312	else if(iflag_ajs.eq.1) then
1313
1314	c Ajustement sec
1315	c ==============
1316	IF(prt_level>9)WRITE(lunout,*)'ajsec'
1317
1318	! ADAPTATION GCM POUR CP(T)
1319	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
1320	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
1321
1322	! ADAPTATION GCM POUR CP(T)
1323	do i=1,klon
1324	flux_ajs(i,1) = 0.0
1325	do j=2,klev
1326	flux_ajs(i,j) = flux_ajs(i,j-1)
1327	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
1328	. *(paprs(i,j-1)-paprs(i,j))
1329	enddo
1330	enddo
1331
1332	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1333	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1334	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1335	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1336	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1337	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1338
1339	if (iflag_trac.eq.1) then
1340	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1341	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1342	endif
1343	endif
1344
1345	! tests: output tendencies
1346	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
1347	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
1348	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
1349	c
1350	IF (if_ebil.ge.2) THEN
1351	ztit='after dry_adjust'
1352	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1353	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1354	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1355	call diagphy(cell_area,ztit,ip_ebil
1356	e , zero_v, zero_v, zero_v, zero_v, sens
1357	e , zero_v, zero_v, zero_v, ztsol
1358	e , d_h_vcol, d_qt, d_ec
1359	s , fs_bound, fq_bound )
1360	END IF
1361
1362	c------------------------------------
1363	c Compute mean mass, cp and R :
1364	c------------------------------------
1365
1366	if(callthermos) then
1367	call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax)
1368	endif
1369
1370	c====================================================================
1371	c RAYONNEMENT
1372	c====================================================================
1373	if (mod(itap,radpas) == 0) then
1374
1375	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1376	c====================================================================
1377
1378	if (callnlte .or. callthermos) then
1379	if (ok_chem) then
1380
1381	! nlte : use computed chemical species
1382
1383	co2vmr_gcm(:,:) = tr_seri(:,:,i_co2)*mmean(:,:)/m_tr(i_co2)
1384	covmr_gcm(:,:) = tr_seri(:,:,i_co)*mmean(:,:)/m_tr(i_co)
1385	ovmr_gcm(:,:) = tr_seri(:,:,i_o)*mmean(:,:)/m_tr(i_o)
1386	n2vmr_gcm(:,:) = tr_seri(:,:,i_n2)*mmean(:,:)/m_tr(i_n2)
1387
1388	else
1389
1390	! nlte : use hedin climatology
1391
1392	call compo_hedin83_mod(pplay,rmu0,
1393	$ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1394	end if
1395	end if
1396
1397	c NLTE cooling from CO2 emission
1398	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1399
1400	IF(callnlte) THEN
1401	if(nltemodel.eq.0.or.nltemodel.eq.1) then
1402	CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
1403	$ tr_seri, d_t_nlte)
1404	else if(nltemodel.eq.2) then
1405	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
1406	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1407	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
1408	if(ierr_nlte.gt.0) then
1409	write(,)
1410	$ 'WARNING: nlte_tcool output with error message',
1411	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
1412	write(,)'I will continue anyway'
1413	endif
1414	endif
1415
1416	ELSE
1417
1418	d_t_nlte(:,:)=0.
1419
1420	ENDIF
1421
1422	c Find number of layers for LTE radiation calculations
1423
1424	IF(callnlte .or. callnirco2)
1425	$ CALL nlthermeq(klon, klev, paprs, pplay)
1426
1427	c
1428	c LTE radiative transfert / solar / IR matrix
1429	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1430	if (physideal) then
1431	CALL radlwsw_newtoncool
1432	e (dist, rmu0, fract, zzlev,
1433	e paprs, pplay,ftsol, t_seri)
1434	else
1435	CALL radlwsw
1436	e (dist, rmu0, fract, zzlev,
1437	e paprs, pplay,ftsol, t_seri)
1438	endif
1439
1440	c albedo variations: test for Yeon Joo Lee
1441	c increment to increase it for 20 Vd => +80%
1442	c heat(:,:)=heat(:,:)(1.+0.80((rjourvrai-356)+gmtime)/20.)
1443	c or to decrease it for 20 Vd => 1/1.8
1444	c heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.)
1445
1446	c CO2 near infrared absorption
1447	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1448
1449	d_t_nirco2(:,:)=0.
1450	if (callnirco2) then
1451	call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
1452	. rmu0, fract, d_t_nirco2)
1453	endif
1454
1455
1456	c Net atmospheric radiative heating rate (K.s-1)
1457	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1458
1459	IF(callnlte.or.callnirco2) THEN
1460	CALL blendrad(klon, klev, pplay,heat,
1461	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
1462	ELSE
1463	dtsw(:,:)=heat(:,:)
1464	dtlw(:,:)=-1*cool(:,:)
1465	ENDIF
1466
1467	DO k=1,klev
1468	DO i=1,klon
1469	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
1470	ENDDO
1471	ENDDO
1472
1473
1474	cc---------------------------------------------
1475
1476	c EUV heating rate (K.s-1)
1477	c ~~~~~~~~~~~~~~~~~~~~~~~~
1478
1479	d_t_euv(:,:)=0.
1480
1481	IF (callthermos) THEN
1482
1483	c call euvheat(klon, klev,t_seri,paprs,pplay,zzlay,
1484	c $ rmu0,pdtphys,gmtime,rjourvrai, co2vmr_gcm, n2vmr_gcm,
1485	c $ covmr_gcm, ovmr_gcm,d_t_euv )
1486	call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
1487	$ rmu0,dtimerad,gmtime,rjourvrai,
1488	$ tr_seri, d_tr, d_t_euv )
1489
1490	DO k=1,klev
1491	DO ig=1,klon
1492	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
1493
1494	ENDDO
1495	ENDDO
1496
1497	ENDIF ! callthermos
1498
1499	c====================================================================
1500	ENDIF ! radpas
1501	c====================================================================
1502	c
1503	c Ajouter la tendance des rayonnements (tous les pas)
1504	c
1505	DO k = 1, klev
1506	DO i = 1, klon
1507	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
1508	ENDDO
1509	ENDDO
1510
1511	! increment physics counter
1512
1513	itap = itap + 1
1514
1515	! CONDUCTION and MOLECULAR VISCOSITY
1516	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1517
1518	d_t_conduc(:,:)=0.
1519	d_u_molvis(:,:)=0.
1520	d_v_molvis(:,:)=0.
1521
1522	IF (callthermos) THEN
1523
1524	tsurf(:)=t_seri(:,1)
1525	call conduction(klon, klev,pdtphys,
1526	$ pplay,paprs,t_seri,
1527	$ tsurf,zzlev,zzlay,d_t_conduc)
1528
1529	call molvis(klon, klev,pdtphys,
1530	$ pplay,paprs,t_seri,
1531	$ u,tsurf,zzlev,zzlay,d_u_molvis)
1532
1533	call molvis(klon, klev, pdtphys,
1534	$ pplay,paprs,t_seri,
1535	$ v,tsurf,zzlev,zzlay,d_v_molvis)
1536
1537	DO k=1,klev
1538	DO ig=1,klon
1539	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
1540	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
1541	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
1542	ENDDO
1543	ENDDO
1544	ENDIF
1545
1546
1547	! -- MOLECULAR DIFFUSION ---
1548
1549	d_q_moldif(:,:,:)=0
1550
1551	IF (callthermos .and. ok_chem) THEN
1552
1553	call moldiff_red(klon, klev, nqmax,
1554	& pplay,paprs,t_seri, tr_seri, pdtphys,
1555	& zzlay,d_t_euv,d_t_conduc,d_q_moldif)
1556
1557
1558	! --- update tendencies tracers ---
1559
1560	DO iq = 1, nqmax
1561	DO k=1,klev
1562	DO ig=1,klon
1563	tr_seri(ig,k,iq)= tr_seri(ig,k,iq)+
1564	& d_q_moldif(ig,k,iq)*dtime ! [Kg/kg]?
1565	ENDDO
1566	ENDDO
1567	ENDDO
1568
1569
1570	ENDIF ! callthermos & ok_chem
1571
1572	c====================================================================
1573	! tests: output tendencies
1574	! call writefield_phy('physiq_dtrad',dtrad,klev)
1575
1576	IF (if_ebil.ge.2) THEN
1577	ztit='after rad'
1578	CALL diagetpq(cell_area,ztit,ip_ebil,2,2,dtime
1579	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1580	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1581	call diagphy(cell_area,ztit,ip_ebil
1582	e , topsw, toplw, solsw, sollw, zero_v
1583	e , zero_v, zero_v, zero_v, ztsol
1584	e , d_h_vcol, d_qt, d_ec
1585	s , fs_bound, fq_bound )
1586	END IF
1587	c
1588
1589	c====================================================================
1590	c Calcul des gravity waves FLOTT
1591	c====================================================================
1592	c
1593	c if (ok_orodr.or.ok_gw_nonoro) then
1594
1595	c CALCUL DE N2
1596	c UTILISE LA RELATION ENTRE N2 ET STABILITE
1597	c N2 = RG/T (dT/dz+RG/cp(T))
1598	c ET DONC EN N'UTILISE QUE LA TEMPERATURE, PAS teta.
1599
1600	do i=1,klon
1601	do k=2,klev
1602	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1603	enddo
1604	enddo
1605	do i=1,klon
1606	do k=2,klev
1607	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1608	zdtlev(i,k) = t_seri(i,k)-t_seri(i,k-1)
1609	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
1610	zn2(i,k) = RG/ztlev(i,k) * ( zdtlev(i,k)/zdzlev(i,k)
1611	. + RG/cpdet(ztlev(i,k)) )
1612	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1613	enddo
1614	zn2(i,1) = 1.e-12 ! securite
1615	enddo
1616
1617	c endif
1618
1619	c ----------------------------ORODRAG
1620	IF (ok_orodr) THEN
1621	c
1622	c selection des points pour lesquels le shema est actif:
1623	igwd=0
1624	DO i=1,klon
1625	itest(i)=0
1626	c IF ((zstd(i).gt.10.0)) THEN
1627	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1628	itest(i)=1
1629	igwd=igwd+1
1630	idx(igwd)=i
1631	ENDIF
1632	ENDDO
1633	c igwdim=MAX(1,igwd)
1634	c
1635	c A ADAPTER POUR VENUS!!! [ TN: c'est fait ! ]
1636	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
1637	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1638	e igwd,idx,itest,
1639	e t_seri, u_seri, v_seri,
1640	s zulow, zvlow, zustrdr, zvstrdr,
1641	s d_t_oro, d_u_oro, d_v_oro,
1642	s zublstrdr,zvblstrdr,znlow,zeff,zbl,
1643	s ztau,tau0,knu2,kbreak)
1644
1645	c print*,"d_u_oro=",d_u_oro(klon/2,:)
1646	c ajout des tendances
1647	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1648	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1649	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1650	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1651	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1652	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1653	c
1654	ELSE
1655	d_t_oro = 0.
1656	d_u_oro = 0.
1657	d_v_oro = 0.
1658	zustrdr = 0.
1659	zvstrdr = 0.
1660	zublstrdr = 0.
1661	zvblstrdr = 0.
1662	znlow = 0.
1663	zeff = 0.
1664	zbl = 0
1665	knu2 = 0
1666	kbreak = 0
1667	ztau = 0
1668	tau0 = 0.
1669	c
1670	ENDIF ! fin de test sur ok_orodr
1671	c
1672	! tests: output tendencies
1673	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
1674	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
1675	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
1676
1677	c ----------------------------OROLIFT
1678	IF (ok_orolf) THEN
1679	print*,"ok_orolf NOT IMPLEMENTED !"
1680	stop
1681	c
1682	c selection des points pour lesquels le shema est actif:
1683	igwd=0
1684	DO i=1,klon
1685	itest(i)=0
1686	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1687	itest(i)=1
1688	igwd=igwd+1
1689	idx(igwd)=i
1690	ENDIF
1691	ENDDO
1692	c igwdim=MAX(1,igwd)
1693	c
1694	c A ADAPTER POUR VENUS!!!
1695	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
1696	c e latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1697	c e igwd,idx,itest,
1698	c e t_seri, u_seri, v_seri,
1699	c s zulow, zvlow, zustrli, zvstrli,
1700	c s d_t_lif, d_u_lif, d_v_lif )
1701
1702	c
1703	c ajout des tendances
1704	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1705	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1706	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1707	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1708	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1709	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1710	c
1711	ELSE
1712	d_t_lif = 0.
1713	d_u_lif = 0.
1714	d_v_lif = 0.
1715	zustrli = 0.
1716	zvstrli = 0.
1717	c
1718	ENDIF ! fin de test sur ok_orolf
1719
1720	c ---------------------------- NON-ORO GRAVITY WAVES
1721	IF(ok_gw_nonoro) then
1722
1723	call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1724	e t_seri, u_seri, v_seri, paprs(klon/2+1,:),
1725	o zustrhi,zvstrhi,
1726	o d_t_hin, d_u_hin, d_v_hin)
1727
1728	c ajout des tendances
1729
1730	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1731	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1732	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1733	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1734	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1735	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1736
1737	ELSE
1738	d_t_hin = 0.
1739	d_u_hin = 0.
1740	d_v_hin = 0.
1741	zustrhi = 0.
1742	zvstrhi = 0.
1743
1744	ENDIF ! fin de test sur ok_gw_nonoro
1745
1746	! tests: output tendencies
1747	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
1748	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
1749	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
1750
1751	c====================================================================
1752	c Transport de ballons
1753	c====================================================================
1754	if (ballons.eq.1) then
1755	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,
1756	& latitude_deg,longitude_deg,
1757	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1758	C t,pplay,u,v,zphi) ! alt above planet average radius
1759	endif !ballons
1760
1761	c====================================================================
1762	c Bilan de mmt angulaire
1763	c====================================================================
1764	if (bilansmc.eq.1) then
1765	CMODDEB FLOTT
1766	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1767	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1768
1769	DO i = 1, klon
1770	zustrph(i)=0.
1771	zvstrph(i)=0.
1772	zustrcl(i)=0.
1773	zvstrcl(i)=0.
1774	ENDDO
1775	DO k = 1, klev
1776	DO i = 1, klon
1777	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1778	c (paprs(i,k)-paprs(i,k+1))/rg
1779	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1780	c (paprs(i,k)-paprs(i,k+1))/rg
1781	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1782	c (paprs(i,k)-paprs(i,k+1))/rg
1783	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1784	c (paprs(i,k)-paprs(i,k+1))/rg
1785	ENDDO
1786	ENDDO
1787
1788	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
1789	C ra,rg,romega,
1790	C latitude_deg,longitude_deg,pphis,
1791	C zustrdr,zustrli,zustrcl,
1792	C zvstrdr,zvstrli,zvstrcl,
1793	C paprs,u,v)
1794
1795	CCMODFIN FLOTT
1796	endif !bilansmc
1797
1798	c====================================================================
1799	c====================================================================
1800	c Calculer le transport de l'eau et de l'energie (diagnostique)
1801	c
1802	c A REVOIR POUR VENUS...
1803	c
1804	c CALL transp (paprs,ftsol,
1805	c e t_seri, q_seri, u_seri, v_seri, zphi,
1806	c s ve, vq, ue, uq)
1807	c
1808	c====================================================================
1809	c+jld ec_conser
1810	DO k = 1, klev
1811	DO i = 1, klon
1812	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1813	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1814	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1815	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1816	END DO
1817	END DO
1818	do i=1,klon
1819	flux_ec(i,1) = 0.0
1820	do j=2,klev
1821	flux_ec(i,j) = flux_ec(i,j-1)
1822	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
1823	enddo
1824	enddo
1825
1826	c-jld ec_conser
1827	c====================================================================
1828	IF (if_ebil.ge.1) THEN
1829	ztit='after physic'
1830	CALL diagetpq(cell_area,ztit,ip_ebil,1,1,dtime
1831	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1832	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1833	C Comme les tendances de la physique sont ajoute dans la dynamique,
1834	C on devrait avoir que la variation d'entalpie par la dynamique
1835	C est egale a la variation de la physique au pas de temps precedent.
1836	C Donc la somme de ces 2 variations devrait etre nulle.
1837	call diagphy(cell_area,ztit,ip_ebil
1838	e , topsw, toplw, solsw, sollw, sens
1839	e , zero_v, zero_v, zero_v, ztsol
1840	e , d_h_vcol, d_qt, d_ec
1841	s , fs_bound, fq_bound )
1842	C
1843	d_h_vcol_phy=d_h_vcol
1844	C
1845	END IF
1846	C
1847	c=======================================================================
1848	c SORTIES
1849	c=======================================================================
1850
1851	c Convertir les incrementations en tendances
1852	c
1853	DO k = 1, klev
1854	DO i = 1, klon
1855	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1856	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1857	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1858	ENDDO
1859	ENDDO
1860	c
1861	DO iq = 1, nqmax
1862	DO k = 1, klev
1863	DO i = 1, klon
1864	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1865	ENDDO
1866	ENDDO
1867	ENDDO
1868
1869	c------------------------
1870	c Calcul moment cinetique
1871	c------------------------
1872	c TEST VENUS...
1873	c mangtot = 0.0
1874	c DO k = 1, klev
1875	c DO i = 1, klon
1876	c mang(i,k) = RA*cos(latitude(i))
1877	c . (u_seri(i,k)+RAcos(latitude(i))*ROMEGA)
1878	c . cell_area(i)(paprs(i,k)-paprs(i,k+1))/RG
1879	c mangtot=mangtot+mang(i,k)
1880	c ENDDO
1881	c ENDDO
1882	c print*,"Moment cinetique total = ",mangtot
1883	c
1884	c------------------------
1885	c
1886	c Sauvegarder les valeurs de t et u a la fin de la physique:
1887	c
1888	DO k = 1, klev
1889	DO i = 1, klon
1890	u_ancien(i,k) = u_seri(i,k)
1891	t_ancien(i,k) = t_seri(i,k)
1892	ENDDO
1893	ENDDO
1894	c
1895	c=============================================================
1896	c Ecriture des sorties
1897	c=============================================================
1898	#ifndef MESOSCALE
1899	#ifdef CPP_IOIPSL
1900
1901	#ifdef histhf
1902	#include "write_histhf.h"
1903	#endif
1904
1905	#ifdef histday
1906	#include "write_histday.h"
1907	#endif
1908
1909	#ifdef histmth
1910	#include "write_histmth.h"
1911	#endif
1912
1913	#ifdef histins
1914	#include "write_histins.h"
1915	#endif
1916
1917	#endif
1918
1919	! XIOS outputs
1920	! This can be done ANYWHERE in the physics routines !
1921
1922	#ifdef CPP_XIOS
1923	! Send fields to XIOS: (NB these fields must also be defined as
1924	! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
1925
1926	! 2D fields
1927
1928	CALL send_xios_field("phis",pphis)
1929	cell_area_out(:)=cell_area(:)
1930	if (is_north_pole_phy) cell_area_out(1)=cell_area(1)/nbp_lon
1931	if (is_south_pole_phy) cell_area_out(klon)=cell_area(klon)/nbp_lon
1932	CALL send_xios_field("aire",cell_area_out)
1933	CALL send_xios_field("tsol",ftsol)
1934	CALL send_xios_field("psol",paprs(:,1))
1935	CALL send_xios_field("cdragh",cdragh)
1936	CALL send_xios_field("cdragm",cdragm)
1937
1938	CALL send_xios_field("tops",topsw)
1939	CALL send_xios_field("topl",toplw)
1940	CALL send_xios_field("sols",solsw)
1941	CALL send_xios_field("soll",sollw)
1942
1943	! 3D fields
1944
1945	CALL send_xios_field("temp",t_seri)
1946	CALL send_xios_field("pres",pplay)
1947	CALL send_xios_field("geop",zphi)
1948	CALL send_xios_field("vitu",u_seri)
1949	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1950	CALL send_xios_field("vitv",-1.*v_seri)
1951	CALL send_xios_field("vitw",omega)
1952	CALL send_xios_field("Kz",ycoefh)
1953	CALL send_xios_field("mmean",mmean)
1954	CALL send_xios_field("rho",rho)
1955	CALL send_xios_field("BV2",zn2)
1956
1957	CALL send_xios_field("dudyn",d_u_dyn)
1958	CALL send_xios_field("duvdf",d_u_vdf)
1959	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1960	CALL send_xios_field("dvvdf",-1.*d_v_vdf)
1961	CALL send_xios_field("duajs",d_u_ajs)
1962	CALL send_xios_field("dugwo",d_u_oro)
1963	CALL send_xios_field("dugwno",d_u_hin)
1964	CALL send_xios_field("dumolvis",d_u_molvis)
1965	c VENUS: regardee a l envers!!!!!!!!!!!!!!!
1966	CALL send_xios_field("dvmolvis",-1.*d_v_molvis)
1967	CALL send_xios_field("dtdyn",d_t_dyn)
1968	CALL send_xios_field("dtphy",d_t)
1969	CALL send_xios_field("dtvdf",d_t_vdf)
1970	CALL send_xios_field("dtajs",d_t_ajs)
1971	CALL send_xios_field("dtswr",dtsw)
1972	CALL send_xios_field("dtswrNLTE",d_t_nirco2)
1973	CALL send_xios_field("dtswrLTE",heat)
1974	CALL send_xios_field("dtlwr",dtlw)
1975	CALL send_xios_field("dtlwrNLTE",d_t_nlte)
1976	CALL send_xios_field("dtlwrLTE",-1.*cool)
1977	CALL send_xios_field("dteuv",d_t_euv)
1978	CALL send_xios_field("dtcond",d_t_conduc)
1979	CALL send_xios_field("dtec",d_t_ec)
1980
1981	CALL send_xios_field("SWnet",swnet(:,1:klev))
1982	CALL send_xios_field("LWnet",lwnet(:,1:klev))
1983	CALL send_xios_field("fluxvdf",fluxt)
1984	CALL send_xios_field("fluxdyn",flux_dyn)
1985	CALL send_xios_field("fluxajs",flux_ajs)
1986	CALL send_xios_field("fluxec",flux_ec)
1987
1988	! when using tracers
1989
1990	if (iflag_trac == 1) then
1991
1992	! tracers in gas phase, volume mixing ratio
1993
1994	do iq = 1,nqmax - nmicro
1995	call send_xios_field(tname(iq),
1996	$ qx(:,:,iq)*mmean(:,:)/m_tr(iq))
1997	end do
1998
1999	! tracers in liquid phase, volume mixing ratio
2000
2001	if ((tr_scheme == 3) .and. (cl_scheme == 1)) THEN
2002	call send_xios_field(tname(i_h2oliq),
2003	$ qx(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
2004	call send_xios_field(tname(i_h2so4liq),
2005	$ qx(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
2006	if (ok_sedim) then
2007	call send_xios_field("Fsedim",fsedim(:,1:klev))
2008	end if
2009	end if
2010
2011	! chemical iterations
2012
2013	call send_xios_field("iter",real(iter))
2014
2015	end if
2016
2017	IF (callthermos .and. ok_chem) THEN
2018	CALL send_xios_field("d_qmoldifCO2",d_q_moldif(:,:,i_co2))
2019	CALL send_xios_field("d_qmoldifO3p",d_q_moldif(:,:,i_o))
2020	CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
2021	ENDIF
2022
2023	if (lafin.and.is_omp_master) then
2024	write(,) "physiq: call xios_context_finalize"
2025	call xios_context_finalize
2026	endif
2027
2028	#endif
2029	#else
2030	! Outputs MESOSCALE
2031	CALL allocate_comm_wrf(klon,klev)
2032	comm_HR_SW(1:klon,1:klev) = dtsw(1:klon,1:klev)
2033	comm_HR_LW(1:klon,1:klev) = dtlw(1:klon,1:klev)
2034	comm_DT_RAD(1:klon,1:klev) = d_t_rad(1:klon,1:klev)
2035	IF (turb_resolved) THEN
2036	open(17,file='hrdyn.txt',form='formatted',status='old')
2037	rewind(17)
2038	DO k=1,klev
2039	read(17,*) dt_dyn(k)
2040	ENDDO
2041	close(17)
2042
2043	do i=1,klon
2044	d_t(i,:)=d_t(i,:)+dt_dyn(:)
2045	comm_HR_DYN(i,:) = dt_dyn(:)
2046	enddo
2047	ELSE
2048	comm_HR_DYN(1:klon,1:klev) = d_t_dyn(1:klon,1:klev)
2049	comm_DT_VDF(1:klon,1:klev) = d_t_vdf(1:klon,1:klev)
2050	comm_DT_AJS(1:klon,1:klev) = d_t_ajs(1:klon,1:klev)
2051	ENDIF
2052	comm_DT(1:klon,1:klev)=d_t(1:klon,1:klev)
2053	#endif
2054
2055
2056	c====================================================================
2057	c Si c'est la fin, il faut conserver l'etat de redemarrage
2058	c====================================================================
2059	c
2060	IF (lafin) THEN
2061	itau_phy = itau_phy + itap
2062	CALL phyredem ("restartphy.nc")
2063
2064	c--------------FLOTT
2065	CMODEB LOTT
2066	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
2067	C DU BILAN DE MOMENT ANGULAIRE.
2068	if (bilansmc.eq.1) then
2069	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
2070	close(27)
2071	close(28)
2072	endif !bilansmc
2073	CMODFIN
2074	c-------------
2075	c--------------SLEBONNOIS
2076	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
2077	C DES BALLONS
2078	if (ballons.eq.1) then
2079	write(,)'Fermeture des ballons*.out'
2080	close(30)
2081	close(31)
2082	close(32)
2083	close(33)
2084	close(34)
2085	endif !ballons
2086	c-------------
2087	ENDIF
2088
2089	END SUBROUTINE physiq
2090
2091	END MODULE physiq_mod
2092

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

Download in other formats:

Original Format