Context Navigation

physiq_mod.F @ 2487

Last change on this file since 2487 was 2486, checked in by emillour, 4 years ago
Venus GCM: Cleanup radlwsw_newtoncool (make it a module in the process) and modify it so that the temperature field towards wich the relmaxation will be done is built using "presnivs" (average model level pressure) and not the GCM's actual pressure field at the first step of the simulation. With this change, one has 1+1=2 when runing with flag "physideal". EM
File size: 67.9 KB

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

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

Download in other formats:

Original Format