Context Navigation

physiq.F @ 1310

Last change on this file since 1310 was 1310, checked in by slebonnois, 10 years ago
SL: VENUS VERTICAL EXTENSION. NLTE and thermospheric processes, to be run with 78 levels and specific inputs.
File size: 50.1 KB

Line
1	!
2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/physiq.F,v 1.8 2005/02/24 09:58:18 fairhead Exp $
3	!
4	c
5	SUBROUTINE physiq (nlon,nlev,nqmax,
6	. debut,lafin,rjourvrai,gmtime,pdtphys,
7	. paprs,pplay,ppk,pphi,pphis,presnivs,
8	. u,v,t,qx,
9	. flxmw,
10	. d_u, d_v, d_t, d_qx, d_ps)
11
12	c======================================================================
13	c
14	c Modifications pour la physique de Venus
15	c S. Lebonnois (LMD/CNRS) Septembre 2005
16	c
17	c ---------------------------------------------------------------------
18	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
19	c
20	c Objet: Moniteur general de la physique du modele
21	cAA Modifications quant aux traceurs :
22	cAA - uniformisation des parametrisations ds phytrac
23	cAA - stockage des moyennes des champs necessaires
24	cAA en mode traceur off-line
25	c modif ( P. Le Van , 12/10/98 )
26	c
27	c Arguments:
28	c
29	c nlon----input-I-nombre de points horizontaux
30	c nlev----input-I-nombre de couches verticales
31	c nqmax---input-I-nombre de traceurs
32	c debut---input-L-variable logique indiquant le premier passage
33	c lafin---input-L-variable logique indiquant le dernier passage
34	c rjour---input-R-numero du jour de l'experience
35	c gmtime--input-R-fraction de la journee (0 a 1)
36	c pdtphys-input-R-pas d'integration pour la physique (seconde)
37	c paprs---input-R-pression pour chaque inter-couche (en Pa)
38	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
39	c ppk ---input-R-fonction d'Exner au milieu de couche
40	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
41	c pphis---input-R-geopotentiel du sol
42	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
43	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
44	c v-------input-R-vitesse Y (de S a N) en m/s
45	c t-------input-R-temperature (K)
46	c qx------input-R-mass mixing ratio traceurs (kg/kg)
47	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
48	c flxmw---input-R-flux de masse vertical en kg/s
49	c
50	c d_u-----output-R-tendance physique de "u" (m/s/s)
51	c d_v-----output-R-tendance physique de "v" (m/s/s)
52	c d_t-----output-R-tendance physique de "t" (K/s)
53	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
54	c d_ps----output-R-tendance physique de la pression au sol
55	c======================================================================
56	USE ioipsl
57	! USE histcom ! not needed; histcom is included in ioipsl
58	USE infotrac
59	USE control_mod
60	use dimphy
61	USE comgeomphy
62	USE mod_phys_lmdz_para, only : is_parallel,jj_nb
63	USE phys_state_var_mod ! Variables sauvegardees de la physique
64	USE write_field_phy
65	USE iophy
66	use cpdet_mod, only: cpdet, t2tpot
67	USE chemparam_mod
68	USE conc
69	USE compo_hedin83_mod2
70	use moyzon_mod, only: tmoy
71	use ieee_arithmetic
72	IMPLICIT none
73	c======================================================================
74	c CLEFS CPP POUR LES IO
75	c =====================
76	c#define histhf
77	#define histday
78	#define histmth
79	#define histins
80	c======================================================================
81	#include "dimensions.h"
82	integer jjmp1
83	parameter (jjmp1=jjm+1-1/jjm)
84	#include "dimsoil.h"
85	#include "clesphys.h"
86	#include "temps.h"
87	#include "iniprint.h"
88	#include "timerad.h"
89	#include "logic.h"
90	#include "tabcontrol.h"
91	#include "nirdata.h"
92	#include "hedin.h"
93	c======================================================================
94	LOGICAL ok_journe ! sortir le fichier journalier
95	save ok_journe
96	c PARAMETER (ok_journe=.true.)
97	c
98	LOGICAL ok_mensuel ! sortir le fichier mensuel
99	save ok_mensuel
100	c PARAMETER (ok_mensuel=.true.)
101	c
102	LOGICAL ok_instan ! sortir le fichier instantane
103	save ok_instan
104	c PARAMETER (ok_instan=.true.)
105	c
106	c======================================================================
107	c
108	c Variables argument:
109	c
110	INTEGER nlon
111	INTEGER nlev
112	INTEGER nqmax
113	REAL rjourvrai
114	REAL gmtime
115	REAL pdtphys
116	LOGICAL debut, lafin
117	REAL paprs(klon,klev+1)
118	REAL pplay(klon,klev)
119	REAL pphi(klon,klev)
120	REAL pphis(klon)
121	REAL presnivs(klev)
122
123	! ADAPTATION GCM POUR CP(T)
124	REAL ppk(klon,klev)
125
126	REAL u(klon,klev)
127	REAL v(klon,klev)
128	REAL t(klon,klev)
129	REAL qx(klon,klev,nqmax)
130
131	REAL d_u_dyn(klon,klev)
132	REAL d_t_dyn(klon,klev)
133
134	REAL flxmw(klon,klev)
135
136	REAL d_u(klon,klev)
137	REAL d_v(klon,klev)
138	REAL d_t(klon,klev)
139	REAL d_qx(klon,klev,nqmax)
140	REAL d_ps(klon)
141
142	logical ok_hf
143	real ecrit_hf
144	integer nid_hf
145	save ok_hf, ecrit_hf, nid_hf
146
147	#ifdef histhf
148	data ok_hf,ecrit_hf/.true.,0.25/
149	#else
150	data ok_hf/.false./
151	#endif
152
153	c Variables propres a la physique
154	c
155	INTEGER,save :: itap ! compteur pour la physique
156	REAL delp(klon,klev) ! epaisseur d'une couche
157	REAL omega(klon,klev) ! vitesse verticale en Pa/s
158
159
160	INTEGER igwd,idx(klon),itest(klon)
161	c
162	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
163
164	REAL zulow(klon),zvlow(klon)
165	REAL zustrdr(klon), zvstrdr(klon)
166	REAL zustrli(klon), zvstrli(klon)
167	REAL zustrhi(klon), zvstrhi(klon)
168
169	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
170	real zdzlev(klon,klev)
171	real ztlev(klon,klev),zpklev(klon,klev)
172	real ztetalay(klon,klev),ztetalev(klon,klev)
173	real zdtetalev(klon,klev)
174	real zn2(klon,klev) ! BV^2 at plev
175
176	c Pour les bilans de moment angulaire,
177	integer bilansmc
178	c Pour le transport de ballons
179	integer ballons
180	c j'ai aussi besoin
181	c du stress de couche limite a la surface:
182
183	REAL zustrcl(klon),zvstrcl(klon)
184
185	c et du stress total c de la physique:
186
187	REAL zustrph(klon),zvstrph(klon)
188
189	c Variables locales:
190	c
191	REAL cdragh(klon) ! drag coefficient pour T and Q
192	REAL cdragm(klon) ! drag coefficient pour vent
193	c
194	cAA Pour TRACEURS
195	cAA
196	REAL,save,allocatable :: source(:,:)
197	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
198	REAL yu1(klon) ! vents dans la premiere couche U
199	REAL yv1(klon) ! vents dans la premiere couche V
200
201	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
202	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
203	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
204	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
205	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
206	c
207	REAL Fsedim(klon,klev+1) ! Flux de sedimentation (kg.m-2)
208
209	c======================================================================
210	c
211	c Declaration des procedures appelees
212	c
213	EXTERNAL ajsec ! ajustement sec
214	EXTERNAL clmain ! couche limite
215	EXTERNAL hgardfou ! verifier les temperatures
216	c EXTERNAL orbite ! calculer l'orbite
217	EXTERNAL phyetat0 ! lire l'etat initial de la physique
218	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
219	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
220	EXTERNAL suphec ! initialiser certaines constantes
221	EXTERNAL transp ! transport total de l'eau et de l'energie
222	EXTERNAL abort_gcm
223	EXTERNAL printflag
224	EXTERNAL zenang
225	EXTERNAL diagetpq
226	EXTERNAL conf_phys
227	EXTERNAL diagphy
228	EXTERNAL mucorr
229	EXTERNAL phytrac
230	EXTERNAL nirco2abs
231	EXTERNAL nir_leedat
232	EXTERNAL nltecool
233	EXTERNAL nlte_tcool
234	EXTERNAL nlte_setup
235	EXTERNAL blendrad
236	EXTERNAL nlthermeq
237	EXTERNAL euvheat
238	EXTERNAL param_read
239	EXTERNAL param_read_e107
240	EXTERNAL conduction
241	EXTERNAL molvis
242	c
243	c Variables locales
244	c
245	CXXX PB
246	REAL fluxt(klon,klev) ! flux turbulent de chaleur
247	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
248	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
249	c
250	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
251	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
252	REAL flux_ec(klon,klev) ! flux de chaleur Ec
253	c
254	REAL tmpout(klon,klev) ! K s-1
255
256	INTEGER itaprad
257	SAVE itaprad
258	c
259	REAL dist, rmu0(klon), fract(klon)
260	REAL zdtime, zlongi
261	c
262	INTEGER i, k, iq, ig, j, ll
263	c
264	REAL zphi(klon,klev)
265	REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
266	real tsurf(klon)
267
268	c va avec nlte_tcool
269	INTEGER ierr_nlte
270	REAL varerr
271
272	c Variables du changement
273	c
274	c ajs: ajustement sec
275	c vdf: couche limite (Vertical DiFfusion)
276	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
277	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
278	c
279	REAL d_ts(klon)
280	c
281	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
282	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
283	c
284	CMOD LOTT: Tendances Orography Sous-maille
285	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
286	REAL d_t_oro(klon,klev)
287	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
288	REAL d_t_lif(klon,klev)
289	C Tendances Ondes de G non oro (runs strato).
290	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
291	REAL d_t_hin(klon,klev)
292
293	c Tendencies due to radiative scheme [K/s]
294	c d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
295	c are not computed at each physical timestep
296	c therefore, they are defined and saved in phys_state_var_mod
297
298	c Tendencies due to molecular viscosity and conduction
299	real d_t_conduc(klon,klev) ! [K/s]
300	real d_u_molvis(klon,klev) ! (m/s) /s
301	real d_v_molvis(klon,klev) ! (m/s) /s
302
303	c
304	c Variables liees a l'ecriture de la bande histoire physique
305	c
306	INTEGER ecrit_mth
307	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
308	c
309	INTEGER ecrit_day
310	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
311	c
312	INTEGER ecrit_ins
313	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
314	c
315	integer itau_w ! pas de temps ecriture = itap + itau_phy
316
317	c Variables locales pour effectuer les appels en serie
318	c
319	REAL t_seri(klon,klev)
320	REAL u_seri(klon,klev), v_seri(klon,klev)
321	c
322	REAL :: tr_seri(klon,klev,nqmax)
323	REAL :: d_tr(klon,klev,nqmax)
324
325	c Variables tendance sedimentation
326
327	REAL :: d_tr_sed(klon,klev,2)
328	REAL :: d_tr_ssed(klon)
329	c
330	c pour ioipsl
331	INTEGER nid_day, nid_mth, nid_ins
332	SAVE nid_day, nid_mth, nid_ins
333	INTEGER nhori, nvert, idayref
334	REAL zsto, zout, zsto1, zsto2, zero
335	parameter (zero=0.0e0)
336	real zjulian
337	save zjulian
338
339	CHARACTER*2 str2
340	character*20 modname
341	character*80 abort_message
342	logical ok_sync
343
344	character*30 nom_fichier
345	character*10 varname
346	character*40 vartitle
347	character*20 varunits
348	C Variables liees au bilan d'energie et d'enthalpi
349	REAL ztsol(klon)
350	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
351	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
352	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
353	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
354	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
355	REAL d_h_vcol_phy
356	REAL fs_bound, fq_bound
357	SAVE d_h_vcol_phy
358	REAL zero_v(klon),zero_v2(klon,klev)
359	CHARACTER*15 ztit
360	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
361	SAVE ip_ebil
362	DATA ip_ebil/2/
363	INTEGER if_ebil ! level for energy conserv. dignostics
364	SAVE if_ebil
365	c+jld ec_conser
366	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
367	c-jld ec_conser
368
369	c TEST VENUS...
370	REAL mang(klon,klev) ! moment cinetique
371	REAL mangtot ! moment cinetique total
372
373	c Declaration des constantes et des fonctions thermodynamiques
374	c
375	#include "YOMCST.h"
376
377	c======================================================================
378	c INITIALISATIONS
379	c================
380
381	modname = 'physiq'
382	ok_sync=.TRUE.
383
384	bilansmc = 0
385	ballons = 0
386	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
387	if (is_parallel) then
388	bilansmc = 0
389	ballons = 0
390	endif
391
392	IF (if_ebil.ge.1) THEN
393	DO i=1,klon
394	zero_v(i)=0.
395	END DO
396	DO i=1,klon
397	DO j=1,klev
398	zero_v2(i,j)=0.
399	END DO
400	END DO
401	END IF
402
403	c PREMIER APPEL SEULEMENT
404	c========================
405	IF (debut) THEN
406	allocate(source(klon,nqmax))
407
408	CALL suphec ! initialiser constantes et parametres phys.
409
410	IF (if_ebil.ge.1) d_h_vcol_phy=0.
411	c
412	c appel a la lecture du physiq.def
413	c
414	call conf_phys(ok_journe, ok_mensuel,
415	. ok_instan,
416	. if_ebil)
417
418	call phys_state_var_init
419	c
420	c Initialising Hedin model for upper atm
421	c (to be revised when coupled to chemistry) :
422	call conc_init
423	c
424	c Initialiser les compteurs:
425	c
426	itap = 0
427	itaprad = 0
428	c
429	c Lecture startphy.nc :
430	c
431	CALL phyetat0 ("startphy.nc")
432
433	c dtime est defini dans tabcontrol.h et lu dans startphy
434	c pdtphys est calcule a partir des nouvelles conditions:
435	c Reinitialisation du pas de temps physique quand changement
436	IF (ABS(dtime-pdtphys).GT.0.001) THEN
437	WRITE(lunout,*) 'Pas physique a change',dtime,
438	. pdtphys
439	c abort_message='Pas physique n est pas correct '
440	c call abort_gcm(modname,abort_message,1)
441	c----------------
442	c pour initialiser convenablement le time_counter, il faut tenir compte
443	c du changement de dtime en changeant itau_phy (point de depart)
444	itau_phy = NINT(itau_phy*dtime/pdtphys)
445	c----------------
446	dtime=pdtphys
447	ENDIF
448
449	radpas = NINT( RDAY/pdtphys/nbapp_rad)
450
451	CALL printflag( ok_journe,ok_instan )
452	c
453	c---------
454	c FLOTT
455	IF (ok_orodr) THEN
456	DO i=1,klon
457	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
458	ENDDO
459	CALL SUGWD(klon,klev,paprs,pplay)
460	DO i=1,klon
461	zuthe(i)=0.
462	zvthe(i)=0.
463	if(zstd(i).gt.10.)then
464	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
465	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
466	endif
467	ENDDO
468	ENDIF
469
470	if (bilansmc.eq.1) then
471	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
472	C DU BILAN DE MOMENT ANGULAIRE.
473	open(27,file='aaam_bud.out',form='formatted')
474	open(28,file='fields_2d.out',form='formatted')
475	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
476	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
477	endif !bilansmc
478
479	c--------------SLEBONNOIS
480	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
481	C DES BALLONS
482	if (ballons.eq.1) then
483	open(30,file='ballons-lat.out',form='formatted')
484	open(31,file='ballons-lon.out',form='formatted')
485	open(32,file='ballons-u.out',form='formatted')
486	open(33,file='ballons-v.out',form='formatted')
487	open(34,file='ballons-alt.out',form='formatted')
488	write(,)'Ouverture des ballons*.out'
489	endif !ballons
490	c-------------
491
492	c---------
493	C TRACEURS
494	C source dans couche limite
495	source = 0.0 ! pas de source, pour l'instant
496	c---------
497
498	c---------
499	c INITIALIZE THERMOSPHERIC PARAMETERS
500	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
501
502	if (callthermos) then
503	if(solvarmod.eq.0) call param_read
504	if(solvarmod.eq.1) call param_read_e107
505	endif
506
507	c Initialisation (recomputed in concentration2)
508	do ig=1,klon
509	do j=1,klev
510	rnew(ig,j)=R
511	cpnew(ig,j)=cpdet(tmoy(j))
512	c print*, ' physique l503'
513	c print*, j, cpdet(tmoy(j))
514	mmean(ig,j)=RMD
515	akknew(ig,j)=1.e-4
516	enddo
517	c stop
518
519	enddo
520
521	IF(callthermos.or.callnlte.or.callnirco2) THEN
522	call compo_hedin83_init2
523	ENDIF
524	if (callnlte) call nlte_setup
525	if(callnirco2.and.(nircorr.eq.1)) call nir_leedat
526	c---------
527
528	c
529	c Verifications:
530	c
531	IF (nlon .NE. klon) THEN
532	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
533	. klon
534	abort_message='nlon et klon ne sont pas coherents'
535	call abort_gcm(modname,abort_message,1)
536	ENDIF
537	IF (nlev .NE. klev) THEN
538	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
539	. klev
540	abort_message='nlev et klev ne sont pas coherents'
541	call abort_gcm(modname,abort_message,1)
542	ENDIF
543	c
544	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
545	$ THEN
546	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
547	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
548	abort_message='Nbre d appels au rayonnement insuffisant'
549	call abort_gcm(modname,abort_message,1)
550	ENDIF
551	c
552	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
553	. iflag_ajs
554	c
555	ecrit_mth = NINT(RDAY/dtime*ecriphy) ! tous les ecritphy jours
556	IF (ok_mensuel) THEN
557	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
558	. ecrit_mth
559	ENDIF
560
561	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
562	IF (ok_journe) THEN
563	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
564	. ecrit_day
565	ENDIF
566
567	ecrit_ins = NINT(RDAY/dtime*ecriphy) ! Fraction de jour reglable
568	IF (ok_instan) THEN
569	WRITE(lunout,*)'La frequence de sortie instant. est de ',
570	. ecrit_ins
571	ENDIF
572
573	c Initialisation des sorties
574	c===========================
575
576	#ifdef CPP_IOIPSL
577
578	#ifdef histhf
579	#include "ini_histhf.h"
580	#endif
581
582	#ifdef histday
583	#include "ini_histday.h"
584	#endif
585
586	#ifdef histmth
587	#include "ini_histmth.h"
588	#endif
589
590	#ifdef histins
591	#include "ini_histins.h"
592	#endif
593
594	#endif
595
596	c
597	c Initialiser les valeurs de u pour calculs tendances
598	c (pour T, c'est fait dans phyetat0)
599	c
600	DO k = 1, klev
601	DO i = 1, klon
602	u_ancien(i,k) = u(i,k)
603	ENDDO
604	ENDDO
605
606	c---------
607	c Ecriture fichier initialisation
608	c PRINT*,'Ecriture Initial_State.csv'
609	c OPEN(88,file='Trac_Point.csv',
610	c & form='formatted')
611	c---------
612
613	c---------
614	c Initialisation des parametres des nuages
615	c===============================================
616
617	if ((nlon .EQ. 1) .AND. ok_cloud) then
618	PRINT*,'Open profile_cloud_parameters.csv'
619	OPEN(66,file='profile_cloud_parameters.csv',
620	& form='formatted')
621	endif
622
623	if ((nlon .EQ. 1) .AND. ok_sedim) then
624	PRINT*,'Open profile_cloud_sedim.csv'
625	OPEN(77,file='profile_cloud_sedim.csv',
626	& form='formatted')
627	endif
628
629	if ((nlon .GT. 1) .AND. ok_chem) then
630	c !!! DONC 3D !!!
631	CALL chemparam_ini()
632	endif
633
634	if ((nlon .GT. 1) .AND. ok_cloud) then
635	c !!! DONC 3D !!!
636	CALL cloud_ini(nlon,nlev)
637	endif
638
639	ENDIF ! debut
640	c====================================================================
641	c======================================================================
642
643	c Mettre a zero des variables de sortie (pour securite)
644	c
645	DO i = 1, klon
646	d_ps(i) = 0.0
647	ENDDO
648	DO k = 1, klev
649	DO i = 1, klon
650	d_t(i,k) = 0.0
651	d_u(i,k) = 0.0
652	d_v(i,k) = 0.0
653	ENDDO
654	ENDDO
655	DO iq = 1, nqmax
656	DO k = 1, klev
657	DO i = 1, klon
658	d_qx(i,k,iq) = 0.0
659	ENDDO
660	ENDDO
661	ENDDO
662	c
663	c Ne pas affecter les valeurs entrees de u, v, h, et q
664	c
665	DO k = 1, klev
666	DO i = 1, klon
667	t_seri(i,k) = t(i,k)
668	u_seri(i,k) = u(i,k)
669	v_seri(i,k) = v(i,k)
670	ENDDO
671	ENDDO
672	DO iq = 1, nqmax
673	DO k = 1, klev
674	DO i = 1, klon
675	tr_seri(i,k,iq) = qx(i,k,iq)
676	ENDDO
677	ENDDO
678	ENDDO
679	C
680	DO i = 1, klon
681	ztsol(i) = ftsol(i)
682	ENDDO
683	C
684	IF (if_ebil.ge.1) THEN
685	ztit='after dynamic'
686	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
687	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
688	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
689	C Comme les tendances de la physique sont ajoute dans la dynamique,
690	C on devrait avoir que la variation d'entalpie par la dynamique
691	C est egale a la variation de la physique au pas de temps precedent.
692	C Donc la somme de ces 2 variations devrait etre nulle.
693	call diagphy(airephy,ztit,ip_ebil
694	e , zero_v, zero_v, zero_v, zero_v, zero_v
695	e , zero_v, zero_v, zero_v, ztsol
696	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
697	s , fs_bound, fq_bound )
698	END IF
699
700	c====================================================================
701	c Diagnostiquer la tendance dynamique
702	c
703	IF (ancien_ok) THEN
704	DO k = 1, klev
705	DO i = 1, klon
706	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
707	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
708	ENDDO
709	ENDDO
710
711	! ADAPTATION GCM POUR CP(T)
712	do i=1,klon
713	flux_dyn(i,1) = 0.0
714	do j=2,klev
715	flux_dyn(i,j) = flux_dyn(i,j-1)
716	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
717	enddo
718	enddo
719
720	ELSE
721	DO k = 1, klev
722	DO i = 1, klon
723	d_u_dyn(i,k) = 0.0
724	d_t_dyn(i,k) = 0.0
725	ENDDO
726	ENDDO
727	ancien_ok = .TRUE.
728	ENDIF
729	c====================================================================
730
731	c Calcule de vitesse verticale a partir de flux de masse verticale
732	DO k = 1, klev
733	DO i = 1, klon
734	omega(i,k) = RG*flxmw(i,k) / airephy(i)
735	END DO
736	END DO
737
738	c
739	c Ajouter le geopotentiel du sol:
740	c
741	DO k = 1, klev
742	DO i = 1, klon
743	zphi(i,k) = pphi(i,k) + pphis(i)
744	ENDDO
745	ENDDO
746
747	c calcul du geopotentiel aux niveaux intercouches
748	c ponderation des altitudes au niveau des couches en dp/p
749
750	DO k=1,klev
751	DO i=1,klon
752	zzlay(i,k)=zphi(i,k)/RG ! [m]
753	ENDDO
754	ENDDO
755	DO i=1,klon
756	zzlev(i,1)=pphis(i)/RG ! [m]
757	ENDDO
758	DO k=2,klev
759	DO i=1,klon
760	z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
761	z2=(paprs(i,k) +pplay(i,k))/(paprs(i,k) -pplay(i,k))
762	zzlev(i,k)=(z1zzlay(i,k-1)+z2zzlay(i,k))/(z1+z2)
763	ENDDO
764	ENDDO
765	DO i=1,klon
766	zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
767	ENDDO
768
769	c====================================================================
770	c
771	c Verifier les temperatures
772	c
773	CALL hgardfou(t_seri,ftsol,'debutphy')
774	c====================================================================
775	c
776	c Incrementer le compteur de la physique
777	c
778	itap = itap + 1
779
780	c====================================================================
781	c Orbite et eclairement
782	c====================================================================
783
784	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
785	c donc pas de variations de Ls, ni de dist.
786	c La seule chose qui compte, c'est la rotation de la planete devant
787	c le Soleil...
788
789	zlongi = 0.0
790	dist = 0.72 ! en UA
791
792	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
793	c a sa valeur, et prendre en compte leur evolution,
794	c il faudra refaire un orbite.F...
795	c CALL orbite(zlongi,dist)
796
797	IF (cycle_diurne) THEN
798	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
799	CALL zenang(zlongi,gmtime,zdtime,rlatd,rlond,rmu0,fract)
800	ELSE
801	call mucorr(klon,zlongi,rlatd,rmu0,fract)
802	ENDIF
803
804	c====================================================================
805	c Calcul des tendances traceurs
806	c====================================================================
807
808	if (iflag_trac.eq.1) then
809
810	if (tr_scheme.eq.1) then
811	! Case 1: pseudo-chemistry with relaxation toward fixed profile
812
813	call phytrac_relax (debut,lafin,nqmax,
814	I nlon,nlev,dtime,pplay,
815	O tr_seri)
816
817	elseif (tr_scheme.eq.2) then
818	! Case 2: surface emission
819	! For the moment, inspired from Mars version
820	! However, the variable 'source' could be used in physiq
821	! so the call to phytrac_emiss could be to initialise it.
822
823	call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
824	I debut,lafin,nqmax,
825	I nlon,nlev,dtime,paprs,
826	I rlatd,rlond,
827	O tr_seri)
828
829	elseif (tr_scheme.eq.3) then ! identical to ok_chem.or.ok_cloud
830	! Case 3: Full chemistry and/or clouds
831
832	call phytrac_chimie(
833	I debut,
834	I gmtime,
835	I nqmax,
836	I nlon,
837	I rlatd,
838	I rlond,
839	I nlev,
840	I dtime,
841	I t_seri,pplay,
842	O tr_seri,
843	O NBRTOT,
844	O WH2SO4,
845	O rho_droplet)
846
847	c CALL WriteField_phy('Pression',pplay,nlev)
848	c CALL WriteField_phy('PressionBnd',paprs,nlev+1)
849	c CALL WriteField_phy('Temp',t_seri,nlev)
850	c IF (ok_cloud) THEN
851	c CALL WriteField_phy('NBRTOT',NBRTOT,nlev)
852	c ENDIF
853	c CALL WriteField_phy('SAl',tr_seri(:,:,i_h2so4liq),nlev)
854	c CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev)
855
856	if (ok_sedim) then
857
858	CALL new_cloud_sedim(
859	I klon,
860	I nlev,
861	I dtime,
862	I pplay,
863	I paprs,
864	I t_seri,
865	I WH2SO4,
866	I tr_seri,
867	I nqmax,
868	I NBRTOT,
869	I rho_droplet,
870	O Fsedim,
871	O d_tr_sed,
872	O d_tr_ssed)
873
874	DO k = 1, klev
875	DO i = 1, klon
876
877	c WRITE(88,"(11(e15.8,','))") pplay(5,25),
878	c & t_seri(5,25),tr_seri(5,25,i_h2oliq),
879	c & tr_seri(5,25,i_h2o),tr_seri(5,25,i_h2so4liq),
880	c & tr_seri(5,25,i_h2so4),NBRTOT(5,25),WH2SO4(5,25),
881	c & Fsedim(5,25),d_tr_sed(5,25,1),d_tr_sed(5,25,2)
882
883	c--------------------
884	c Ce test est necessaire pour eviter Xliq=NaN
885	IF (ieee_is_nan(d_tr_sed(i,k,1)).OR.
886	& ieee_is_nan(d_tr_sed(i,k,2))) THEN
887	PRINT*,'sedim NaN PROBLEM'
888	PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k)
889	PRINT*,'lat-lon',i,'level',k,'dtime',dtime
890	PRINT*,'NBRTOT',NBRTOT(i,k),'F_sed',Fsedim(i,k)
891	PRINT*,'==============================================='
892	d_tr_sed(i,k,:)=0.
893	ENDIF
894	c--------------------
895
896	tr_seri(i,k,i_h2so4liq) = tr_seri(i,k,i_h2so4liq)+
897	& d_tr_sed(i,k,1)
898	tr_seri(i,k,i_h2oliq) = tr_seri(i,k,i_h2oliq)+
899	& d_tr_sed(i,k,2)
900	d_tr_sed(i,k,:) = d_tr_sed(i,k,:) / dtime
901	Fsedim(i,k) = Fsedim(i,k) / dtime
902
903	ENDDO
904	ENDDO
905
906	Fsedim(:,klev+1) = 0.
907
908	endif ! ok_sedim
909
910	endif ! tr_scheme
911	endif ! iflag_trac
912
913	c
914	c====================================================================
915	c Appeler la diffusion verticale (programme de couche limite)
916	c====================================================================
917
918	c-------------------------------
919	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
920	c l'equilibre radiatif du sol
921	if (1.eq.0) then
922	if (debut) then
923	print*,"ATTENTION, CLMAIN SHUNTEE..."
924	endif
925
926	DO i = 1, klon
927	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
928	fder(i) = 0.0e0
929	dlw(i) = 0.0e0
930	ENDDO
931
932	c Incrementer la temperature du sol
933	c
934	DO i = 1, klon
935	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
936	ftsol(i) = ftsol(i) + d_ts(i)
937	do j=1,nsoilmx
938	ftsoil(i,j)=ftsol(i)
939	enddo
940	ENDDO
941
942	c-------------------------------
943	else
944	c-------------------------------
945
946	fder = dlw
947
948	! ADAPTATION GCM POUR CP(T)
949
950	CALL clmain(dtime,itap,
951	e t_seri,u_seri,v_seri,
952	e rmu0,
953	e ftsol,
954	$ ftsoil,
955	$ paprs,pplay,ppk,radsol,falbe,
956	e solsw, sollw, sollwdown, fder,
957	e rlond, rlatd, cuphy, cvphy,
958	e debut, lafin,
959	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
960	s fluxt,fluxu,fluxv,cdragh,cdragm,
961	s dsens,
962	s ycoefh,yu1,yv1)
963
964	CXXX Incrementation des flux
965	DO i = 1, klon
966	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
967	fder(i) = dlw(i) + dsens(i)
968	ENDDO
969	CXXX
970
971	DO k = 1, klev
972	DO i = 1, klon
973	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
974	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
975	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
976	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
977	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
978	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
979	ENDDO
980	ENDDO
981
982	C TRACEURS
983
984	if (iflag_trac.eq.1) then
985	DO k = 1, klev
986	DO i = 1, klon
987	delp(i,k) = paprs(i,k)-paprs(i,k+1)
988	ENDDO
989	ENDDO
990	DO iq=1, nqmax
991	CALL cltrac(dtime,ycoefh,t_seri,
992	s tr_seri(1,1,iq),source,
993	e paprs, pplay,delp,
994	s d_tr_vdf(1,1,iq))
995	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
996	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
997	ENDDO
998	endif
999
1000	IF (if_ebil.ge.2) THEN
1001	ztit='after clmain'
1002	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
1003	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1004	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1005	call diagphy(airephy,ztit,ip_ebil
1006	e , zero_v, zero_v, zero_v, zero_v, sens
1007	e , zero_v, zero_v, zero_v, ztsol
1008	e , d_h_vcol, d_qt, d_ec
1009	s , fs_bound, fq_bound )
1010	END IF
1011	C
1012	c
1013	c Incrementer la temperature du sol
1014	c
1015	DO i = 1, klon
1016	ftsol(i) = ftsol(i) + d_ts(i)
1017	ENDDO
1018
1019	c Calculer la derive du flux infrarouge
1020	c
1021	DO i = 1, klon
1022	dlw(i) = - 4.0RSIGMAftsol(i)**3
1023	ENDDO
1024
1025	c-------------------------------
1026	endif ! fin du VENUS TEST
1027
1028	! tests: output tendencies
1029	! call writefield_phy('physiq_d_t_vdf',d_t_vdf,klev)
1030	! call writefield_phy('physiq_d_u_vdf',d_u_vdf,klev)
1031	! call writefield_phy('physiq_d_v_vdf',d_v_vdf,klev)
1032	! call writefield_phy('physiq_d_ts',d_ts,1)
1033
1034	c
1035	c Appeler l'ajustement sec
1036	c
1037	c===================================================================
1038	c Convection seche
1039	c===================================================================
1040	c
1041	d_t_ajs(:,:)=0.
1042	d_u_ajs(:,:)=0.
1043	d_v_ajs(:,:)=0.
1044	d_tr_ajs(:,:,:)=0.
1045	c
1046	IF(prt_level>9)WRITE(lunout,*)
1047	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
1048	s ,iflag_ajs
1049
1050	if(iflag_ajs.eq.0) then
1051	c Rien
1052	c ====
1053	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1054
1055	else if(iflag_ajs.eq.1) then
1056
1057	c Ajustement sec
1058	c ==============
1059	IF(prt_level>9)WRITE(lunout,*)'ajsec'
1060
1061	! ADAPTATION GCM POUR CP(T)
1062	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
1063	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
1064
1065	! ADAPTATION GCM POUR CP(T)
1066	do i=1,klon
1067	flux_ajs(i,1) = 0.0
1068	do j=2,klev
1069	flux_ajs(i,j) = flux_ajs(i,j-1)
1070	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
1071	. *(paprs(i,j-1)-paprs(i,j))
1072	enddo
1073	enddo
1074
1075	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1076	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
1077	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
1078	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
1079	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
1080	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
1081
1082	if (iflag_trac.eq.1) then
1083	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
1084	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
1085	endif
1086
1087	endif
1088
1089	! tests: output tendencies
1090	! call writefield_phy('physiq_d_t_ajs',d_t_ajs,klev)
1091	! call writefield_phy('physiq_d_u_ajs',d_u_ajs,klev)
1092	! call writefield_phy('physiq_d_v_ajs',d_v_ajs,klev)
1093	c
1094	IF (if_ebil.ge.2) THEN
1095	ztit='after dry_adjust'
1096	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1097	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1098	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1099	call diagphy(airephy,ztit,ip_ebil
1100	e , zero_v, zero_v, zero_v, zero_v, sens
1101	e , zero_v, zero_v, zero_v, ztsol
1102	e , d_h_vcol, d_qt, d_ec
1103	s , fs_bound, fq_bound )
1104	END IF
1105
1106	c====================================================================
1107	c RAYONNEMENT
1108	c====================================================================
1109
1110	c-----------------------------------------------------------------------
1111	c . Compute radiative tendencies :
1112	c------------------------------------
1113	c====================================================================
1114	IF (MOD(itaprad,radpas).EQ.0) THEN
1115	c====================================================================
1116
1117	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1118	c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
1119
1120	c Calcul pour Cp rnew et mmean avec traceurs (Cp independant de T !! )
1121	IF(callnlte.or.callthermos) THEN
1122	call compo_hedin83_mod(pplay,rmu0,
1123	& co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
1124	ENDIF
1125
1126	if(callthermos) then
1127	call concentrations2(pplay,t_seri,d_t,co2vmr_gcm, n2vmr_gcm,
1128	& covmr_gcm, ovmr_gcm,nvmr_gcm,pdtphys)
1129	endif
1130
1131
1132	c NLTE cooling from CO2 emission
1133	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1134
1135	IF(callnlte) THEN
1136	if(nltemodel.eq.0.or.nltemodel.eq.1) then
1137	CALL nltecool(klon, klev, pplay*9.869e-6, t_seri,
1138	$ co2vmr_gcm,n2vmr_gcm, covmr_gcm, ovmr_gcm,
1139	$ d_t_nlte)
1140	else if(nltemodel.eq.2) then
1141	CALL nlte_tcool(klon,klev,pplay*9.869e-6,
1142	$ t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1143	$ ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
1144	if(ierr_nlte.gt.0) then
1145	write(,)
1146	$ 'WARNING: nlte_tcool output with error message',
1147	$ 'ierr_nlte=',ierr_nlte,'varerr=',varerr
1148	write(,)'I will continue anyway'
1149	endif
1150
1151	endif
1152
1153	ELSE
1154
1155	d_t_nlte(:,:)=0.
1156
1157	ENDIF
1158
1159	c Find number of layers for LTE radiation calculations
1160
1161	IF(callnlte .or. callnirco2)
1162	$ CALL nlthermeq(klon, klev, paprs, pplay)
1163
1164
1165	c LTE radiative transfert / solar / IR matrix
1166	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1167
1168	CALL radlwsw
1169	e (dist, rmu0, fract, zzlev,
1170	e paprs, pplay,ftsol, t_seri)
1171
1172
1173	c CO2 near infrared absorption
1174	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1175
1176	d_t_nirco2(:,:)=0.
1177	if (callnirco2) then
1178	call nirco2abs (klon, klev, pplay, dist, nqmax, qx,
1179	. rmu0, fract, d_t_nirco2,
1180	. co2vmr_gcm, ovmr_gcm)
1181	endif
1182
1183
1184	c Net atmospheric radiative heating rate (K.s-1)
1185	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1186
1187	IF(callnlte.or.callnirco2) THEN
1188	CALL blendrad(klon, klev, pplay,heat,
1189	& cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
1190	ELSE
1191	dtsw(:,:)=heat(:,:)
1192	dtlw(:,:)=-1*cool(:,:)
1193	ENDIF
1194
1195	DO k=1,klev
1196	DO i=1,klon
1197	d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k) ! K/s
1198	ENDDO
1199	ENDDO
1200
1201
1202	cc---------------------------------------------
1203
1204	c EUV heating rate (K.s-1)
1205	c ~~~~~~~~~~~~~~~~~~~~~~~~
1206
1207	d_t_euv(:,:)=0.
1208
1209	IF (callthermos) THEN
1210
1211	c call thermosphere(zplev,zplay,dist_sol,
1212	c $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
1213	c & pt,pq,pu,pv,pdt,pdq,
1214	c $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
1215
1216	call euvheat(klon, klev, t_seri,paprs,pplay,zzlay,
1217	$ rmu0,pdtphys,gmtime,rjourvrai,
1218	C $ pq,pdq,zdteuv)
1219	$ co2vmr_gcm, n2vmr_gcm, covmr_gcm,
1220	$ ovmr_gcm,nvmr_gcm,d_t_euv )
1221
1222	DO k=1,klev
1223	DO ig=1,klon
1224	d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
1225
1226	ENDDO
1227	ENDDO
1228
1229	ENDIF ! callthermos
1230
1231	c====================================================================
1232	itaprad = 0
1233	ENDIF ! radpas
1234	c====================================================================
1235	c
1236	c Ajouter la tendance des rayonnements (tous les pas)
1237	c
1238	DO k = 1, klev
1239	DO i = 1, klon
1240	t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
1241	ENDDO
1242	ENDDO
1243
1244	! CONDUCTION and MOLECULAR VISCOSITY
1245	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1246
1247	d_t_conduc(:,:)=0.
1248	d_u_molvis(:,:)=0.
1249	d_v_molvis(:,:)=0.
1250
1251	IF (callthermos) THEN
1252
1253	tsurf(:)=t_seri(:,1)
1254	call conduction(klon, klev,pdtphys,
1255	$ pplay,paprs,t_seri,
1256	$ tsurf,zzlev,zzlay,d_t_conduc)
1257
1258	call molvis(klon, klev,pdtphys,
1259	$ pplay,paprs,t_seri,
1260	$ u,tsurf,zzlev,zzlay,d_u_molvis)
1261
1262	call molvis(klon, klev, pdtphys,
1263	$ pplay,paprs,t_seri,
1264	$ v,tsurf,zzlev,zzlay,d_u_molvis)
1265
1266	DO k=1,klev
1267	DO ig=1,klon
1268	t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
1269	u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
1270	v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
1271	ENDDO
1272	ENDDO
1273
1274	ENDIF ! callthermos
1275
1276	c====================================================================
1277	! tests: output tendencies
1278	! call writefield_phy('physiq_dtrad',dtrad,klev)
1279
1280	IF (if_ebil.ge.2) THEN
1281	ztit='after rad'
1282	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1283	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1284	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1285	call diagphy(airephy,ztit,ip_ebil
1286	e , topsw, toplw, solsw, sollw, zero_v
1287	e , zero_v, zero_v, zero_v, ztsol
1288	e , d_h_vcol, d_qt, d_ec
1289	s , fs_bound, fq_bound )
1290	END IF
1291	c
1292
1293	c====================================================================
1294	c Calcul des gravity waves FLOTT
1295	c====================================================================
1296	c
1297	if (ok_orodr.or.ok_gw_nonoro) then
1298	c CALCUL DE N2
1299	do i=1,klon
1300	do k=2,klev
1301	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
1302	zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
1303	enddo
1304	enddo
1305	call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
1306	call t2tpot(klon*klev,t_seri,ztetalay,ppk)
1307	do i=1,klon
1308	do k=2,klev
1309	zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
1310	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
1311	zn2(i,k) = RGzdtetalev(i,k)/(ztetalev(i,k)zdzlev(i,k))
1312	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
1313	enddo
1314	zn2(i,1) = 1.e-12 ! securite
1315	enddo
1316
1317	endif
1318
1319	c ----------------------------ORODRAG
1320	IF (ok_orodr) THEN
1321	c
1322	c selection des points pour lesquels le shema est actif:
1323	igwd=0
1324	DO i=1,klon
1325	itest(i)=0
1326	c IF ((zstd(i).gt.10.0)) THEN
1327	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1328	itest(i)=1
1329	igwd=igwd+1
1330	idx(igwd)=i
1331	ENDIF
1332	ENDDO
1333	c igwdim=MAX(1,igwd)
1334	c
1335	c A ADAPTER POUR VENUS!!!
1336	CALL drag_noro(klon,klev,dtime,paprs,pplay,pphi,zn2,
1337	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
1338	e igwd,idx,itest,
1339	e t_seri, u_seri, v_seri,
1340	s zulow, zvlow, zustrdr, zvstrdr,
1341	s d_t_oro, d_u_oro, d_v_oro)
1342
1343	c print*,"d_u_oro=",d_u_oro(klon/2,:)
1344	c ajout des tendances
1345	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
1346	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
1347	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
1348	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
1349	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
1350	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
1351	c
1352	ELSE
1353	d_t_oro = 0.
1354	d_u_oro = 0.
1355	d_v_oro = 0.
1356	zustrdr = 0.
1357	zvstrdr = 0.
1358	c
1359	ENDIF ! fin de test sur ok_orodr
1360	c
1361	! tests: output tendencies
1362	! call writefield_phy('physiq_d_t_oro',d_t_oro,klev)
1363	! call writefield_phy('physiq_d_u_oro',d_u_oro,klev)
1364	! call writefield_phy('physiq_d_v_oro',d_v_oro,klev)
1365
1366	c ----------------------------OROLIFT
1367	IF (ok_orolf) THEN
1368	print*,"ok_orolf NOT IMPLEMENTED !"
1369	stop
1370	c
1371	c selection des points pour lesquels le shema est actif:
1372	igwd=0
1373	DO i=1,klon
1374	itest(i)=0
1375	IF ((zpic(i)-zmea(i)).GT.100.) THEN
1376	itest(i)=1
1377	igwd=igwd+1
1378	idx(igwd)=i
1379	ENDIF
1380	ENDDO
1381	c igwdim=MAX(1,igwd)
1382	c
1383	c A ADAPTER POUR VENUS!!!
1384	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
1385	c e rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
1386	c e igwd,idx,itest,
1387	c e t_seri, u_seri, v_seri,
1388	c s zulow, zvlow, zustrli, zvstrli,
1389	c s d_t_lif, d_u_lif, d_v_lif )
1390
1391	c
1392	c ajout des tendances
1393	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
1394	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
1395	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
1396	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
1397	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
1398	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
1399	c
1400	ELSE
1401	d_t_lif = 0.
1402	d_u_lif = 0.
1403	d_v_lif = 0.
1404	zustrli = 0.
1405	zvstrli = 0.
1406	c
1407	ENDIF ! fin de test sur ok_orolf
1408
1409	c ---------------------------- NON-ORO GRAVITY WAVES
1410	IF(ok_gw_nonoro) then
1411
1412	call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
1413	e t_seri, u_seri, v_seri,
1414	o zustrhi,zvstrhi,
1415	o d_t_hin, d_u_hin, d_v_hin)
1416
1417	c ajout des tendances
1418
1419	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
1420	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
1421	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
1422	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
1423	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
1424	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
1425
1426	ELSE
1427	d_t_hin = 0.
1428	d_u_hin = 0.
1429	d_v_hin = 0.
1430	zustrhi = 0.
1431	zvstrhi = 0.
1432
1433	ENDIF ! fin de test sur ok_gw_nonoro
1434
1435	! tests: output tendencies
1436	! call writefield_phy('physiq_d_t_hin',d_t_hin,klev)
1437	! call writefield_phy('physiq_d_u_hin',d_u_hin,klev)
1438	! call writefield_phy('physiq_d_v_hin',d_v_hin,klev)
1439
1440	c====================================================================
1441	c Transport de ballons
1442	c====================================================================
1443	if (ballons.eq.1) then
1444	CALL ballon(30,pdtphys,rjourvrai,gmtime*RDAY,rlatd,rlond,
1445	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
1446	C t,pplay,u,v,zphi) ! alt above planet average radius
1447	endif !ballons
1448
1449	c====================================================================
1450	c Bilan de mmt angulaire
1451	c====================================================================
1452	if (bilansmc.eq.1) then
1453	CMODDEB FLOTT
1454	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
1455	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
1456
1457	DO i = 1, klon
1458	zustrph(i)=0.
1459	zvstrph(i)=0.
1460	zustrcl(i)=0.
1461	zvstrcl(i)=0.
1462	ENDDO
1463	DO k = 1, klev
1464	DO i = 1, klon
1465	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
1466	c (paprs(i,k)-paprs(i,k+1))/rg
1467	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
1468	c (paprs(i,k)-paprs(i,k+1))/rg
1469	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
1470	c (paprs(i,k)-paprs(i,k+1))/rg
1471	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
1472	c (paprs(i,k)-paprs(i,k+1))/rg
1473	ENDDO
1474	ENDDO
1475
1476	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime*RDAY,
1477	C ra,rg,romega,
1478	C rlatd,rlond,pphis,
1479	C zustrdr,zustrli,zustrcl,
1480	C zvstrdr,zvstrli,zvstrcl,
1481	C paprs,u,v)
1482
1483	CCMODFIN FLOTT
1484	endif !bilansmc
1485
1486	c====================================================================
1487	c====================================================================
1488	c Calculer le transport de l'eau et de l'energie (diagnostique)
1489	c
1490	c A REVOIR POUR VENUS...
1491	c
1492	c CALL transp (paprs,ftsol,
1493	c e t_seri, q_seri, u_seri, v_seri, zphi,
1494	c s ve, vq, ue, uq)
1495	c
1496	c====================================================================
1497	c+jld ec_conser
1498	DO k = 1, klev
1499	DO i = 1, klon
1500	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
1501	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
1502	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
1503	d_t_ec(i,k) = d_t_ec(i,k)/dtime
1504	END DO
1505	END DO
1506	do i=1,klon
1507	flux_ec(i,1) = 0.0
1508	do j=2,klev
1509	flux_ec(i,j) = flux_ec(i,j-1)
1510	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
1511	enddo
1512	enddo
1513
1514	c-jld ec_conser
1515	c====================================================================
1516	IF (if_ebil.ge.1) THEN
1517	ztit='after physic'
1518	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1519	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
1520	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1521	C Comme les tendances de la physique sont ajoute dans la dynamique,
1522	C on devrait avoir que la variation d'entalpie par la dynamique
1523	C est egale a la variation de la physique au pas de temps precedent.
1524	C Donc la somme de ces 2 variations devrait etre nulle.
1525	call diagphy(airephy,ztit,ip_ebil
1526	e , topsw, toplw, solsw, sollw, sens
1527	e , zero_v, zero_v, zero_v, ztsol
1528	e , d_h_vcol, d_qt, d_ec
1529	s , fs_bound, fq_bound )
1530	C
1531	d_h_vcol_phy=d_h_vcol
1532	C
1533	END IF
1534	C
1535	c=======================================================================
1536	c SORTIES
1537	c=======================================================================
1538
1539	c Convertir les incrementations en tendances
1540	c
1541	DO k = 1, klev
1542	DO i = 1, klon
1543	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
1544	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
1545	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
1546	ENDDO
1547	ENDDO
1548	c
1549	DO iq = 1, nqmax
1550	DO k = 1, klev
1551	DO i = 1, klon
1552	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
1553	ENDDO
1554	ENDDO
1555	ENDDO
1556
1557	c------------------------
1558	c Calcul moment cinetique
1559	c------------------------
1560	c TEST VENUS...
1561	c mangtot = 0.0
1562	c DO k = 1, klev
1563	c DO i = 1, klon
1564	c mang(i,k) = RAcos(rlatd(i)RPI/180.)
1565	c . (u_seri(i,k)+RAcos(rlatd(i)RPI/180.)ROMEGA)
1566	c . airephy(i)(paprs(i,k)-paprs(i,k+1))/RG
1567	c mangtot=mangtot+mang(i,k)
1568	c ENDDO
1569	c ENDDO
1570	c print*,"Moment cinetique total = ",mangtot
1571	c
1572	c------------------------
1573	c
1574	c Sauvegarder les valeurs de t et u a la fin de la physique:
1575	c
1576	DO k = 1, klev
1577	DO i = 1, klon
1578	u_ancien(i,k) = u_seri(i,k)
1579	t_ancien(i,k) = t_seri(i,k)
1580	ENDDO
1581	ENDDO
1582	c
1583	c=============================================================
1584	c Ecriture des sorties
1585	c=============================================================
1586
1587	#ifdef CPP_IOIPSL
1588
1589	#ifdef histhf
1590	#include "write_histhf.h"
1591	#endif
1592
1593	#ifdef histday
1594	#include "write_histday.h"
1595	#endif
1596
1597	#ifdef histmth
1598	#include "write_histmth.h"
1599	#endif
1600
1601	#ifdef histins
1602	#include "write_histins.h"
1603	#endif
1604
1605	#endif
1606
1607	c====================================================================
1608	c Si c'est la fin, il faut conserver l'etat de redemarrage
1609	c====================================================================
1610	c
1611	IF (lafin) THEN
1612	itau_phy = itau_phy + itap
1613	CALL phyredem ("restartphy.nc")
1614
1615	c--------------FLOTT
1616	CMODEB LOTT
1617	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
1618	C DU BILAN DE MOMENT ANGULAIRE.
1619	if (bilansmc.eq.1) then
1620	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
1621	close(27)
1622	close(28)
1623	endif !bilansmc
1624	CMODFIN
1625	c-------------
1626	c--------------SLEBONNOIS
1627	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
1628	C DES BALLONS
1629	if (ballons.eq.1) then
1630	write(,)'Fermeture des ballons*.out'
1631	close(30)
1632	close(31)
1633	close(32)
1634	close(33)
1635	close(34)
1636	endif !ballons
1637	c-------------
1638	ENDIF
1639
1640	RETURN
1641	END
1642
1643
1644
1645	***********************************************************************
1646	***********************************************************************
1647	***********************************************************************
1648	***********************************************************************
1649	***********************************************************************
1650	***********************************************************************
1651	***********************************************************************
1652	***********************************************************************
1653
1654	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
1655	IMPLICIT none
1656	c
1657	c Tranformer une variable de la grille physique a
1658	c la grille d'ecriture
1659	c
1660	INTEGER nfield,nlon,iim,jjmp1, jjm
1661	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
1662	c
1663	INTEGER i, n, ig
1664	c
1665	jjm = jjmp1 - 1
1666	DO n = 1, nfield
1667	DO i=1,iim
1668	ecrit(i,n) = fi(1,n)
1669	ecrit(i+jjm*iim,n) = fi(nlon,n)
1670	ENDDO
1671	DO ig = 1, nlon - 2
1672	ecrit(iim+ig,n) = fi(1+ig,n)
1673	ENDDO
1674	ENDDO
1675	RETURN
1676	END

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

Download in other formats:

Original Format