Context Navigation

physiq.F @ 1452

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

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

Download in other formats:

Original Format