Context Navigation

leapfrog_nogcm_pluto.F @ 3567

Last change on this file since 3567 was 3229, checked in by tbertrand, 21 months ago
Pluto PCM: Adding Nogcm.F : 2D surface model with 1 atm layer
File size: 32.7 KB

Line
1	!
2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
3	!
4	c
5	c
6	SUBROUTINE leapfrog_nogcm_pluto(ucov,vcov,teta,ps,masse,phis,q,
7	& time_0)
8
9
10	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
11	#ifdef CPP_IOIPSL
12	use IOIPSL
13	#endif
14	USE infotrac, ONLY: nqtot,ok_iso_verif,tname
15	USE guide_mod, ONLY : guide_main
16	USE write_field, ONLY: writefield
17	USE control_mod, ONLY: planet_type,nday,day_step,iperiod,iphysiq,
18	& less1day,fractday,ndynstep,iconser,
19	& dissip_period,offline,ip_ebil_dyn,
20	& ok_dynzon,periodav,ok_dyn_ave,iecri,
21	& ok_dyn_ins,output_grads_dyn
22	use exner_hyb_m, only: exner_hyb
23	use exner_milieu_m, only: exner_milieu
24	use cpdet_mod, only: cpdet,tpot2t,t2tpot
25	use sponge_mod, only: callsponge,mode_sponge,sponge
26	use comuforc_h
27	USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs,
28	& aps,bps,presnivs,pseudoalt,preff,scaleheight
29	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,
30	. cpp,ihf,iflag_top_bound,pi,kappa,r
31	USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
32	. statcl,conser,purmats,tidal,ok_strato
33	USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
34	. start_time,dt
35
36
37
38
39	IMPLICIT NONE
40
41	c ...... Version du 10/01/98 ..........
42
43	c avec coordonnees verticales hybrides
44	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
45
46	c=======================================================================
47	c
48	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
49	c -------
50	c
51	c Objet:
52	c ------
53	c
54	c GCM LMD nouvelle grille
55	c
56	c=======================================================================
57	c
58	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
59	c et possibilite d'appeler une fonction f(y) a derivee tangente
60	c hyperbolique a la place de la fonction a derivee sinusoidale.
61
62	c ... Possibilite de choisir le shema pour l'advection de
63	c q , en modifiant iadv dans traceur.def (10/02) .
64	c
65	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
66	c Pour Van-Leer iadv=10
67	c
68	c-----------------------------------------------------------------------
69	c Declarations:
70	c -------------
71
72	#include "dimensions.h"
73	#include "paramet.h"
74	#include "comdissnew.h"
75	#include "comgeom.h"
76	!#include "com_io_dyn.h"
77	#include "iniprint.h"
78	#include "academic.h"
79
80	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
81	! #include "clesphys.h"
82
83	REAL,INTENT(IN) :: time_0 ! not used
84
85	c dynamical variables:
86	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
87	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
88	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
89	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
90	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
91	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
92	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
93
94	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
95	REAL pks(ip1jmp1) ! exner at the surface
96	REAL pk(ip1jmp1,llm) ! exner at mid-layer
97	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
98	REAL phi(ip1jmp1,llm) ! geopotential
99	REAL w(ip1jmp1,llm) ! vertical velocity
100	REAL kpd(ip1jmp1) ! TB15 = exp (-z/H)
101
102	! ADAPTATION GCM POUR CP(T)
103	REAL temp(ip1jmp1,llm) ! temperature
104	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
105
106	! real zqmin,zqmax
107
108
109	! ED18 nogcm
110	REAL tau_ps
111	REAL tau_n2
112	REAL tau_teta
113	REAL tetadpmean
114	REAL tetadp(ip1jmp1,llm)
115	REAL dpmean
116	REAL dp(ip1jmp1,llm)
117	REAL tetamean
118	real globaverage2d
119	! LOGICAL firstcall_globaverage2d
120
121
122	c variables dynamiques intermediaire pour le transport
123	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
124
125	c variables dynamiques au pas -1
126	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
127	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
128	REAL massem1(ip1jmp1,llm)
129
130	c tendances dynamiques en */s
131	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
132	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot)
133
134	c tendances de la dissipation en */s
135	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
136	REAL dtetadis(ip1jmp1,llm)
137
138	c tendances de la couche superieure */s
139	c REAL dvtop(ip1jm,llm)
140	REAL dutop(ip1jmp1,llm)
141	c REAL dtetatop(ip1jmp1,llm)
142	c REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
143
144	c TITAN : tendances due au forces de marees */s
145	REAL dvtidal(ip1jm,llm),dutidal(ip1jmp1,llm)
146
147	c tendances physiques */s
148	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
149	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
150
151	c variables pour le fichier histoire
152	! REAL dtav ! intervalle de temps elementaire
153
154	REAL tppn(iim),tpps(iim),tpn,tps
155	c
156	INTEGER itau,itaufinp1,iav
157	! INTEGER iday ! jour julien
158	REAL time
159
160	REAL SSUM
161	! REAL finvmaold(ip1jmp1,llm)
162
163	cym LOGICAL lafin
164	LOGICAL :: lafin=.false.
165	INTEGER ij,iq,l
166	! INTEGER ik
167
168	! real time_step, t_wrt, t_ops
169
170	REAL rdaym_ini
171	! jD_cur: jour julien courant
172	! jH_cur: heure julienne courante
173	REAL :: jD_cur, jH_cur
174	! INTEGER :: an, mois, jour
175	! REAL :: secondes
176
177	LOGICAL first,callinigrads
178	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
179	save first
180	data first/.true./
181	! real dt_cum
182	! character*10 infile
183	! integer zan, tau0, thoriid
184	! integer nid_ctesGCM
185	! save nid_ctesGCM
186	! real degres
187	! real rlong(iip1), rlatg(jjp1)
188	! real zx_tmp_2d(iip1,jjp1)
189	! integer ndex2d(iip1*jjp1)
190	logical ok_sync
191	parameter (ok_sync = .true.)
192	logical physics
193
194	data callinigrads/.true./
195	character*10 string10
196
197	! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
198	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
199
200	REAL psmean ! pression moyenne
201	REAL pqn2mean ! moyenne globale ps*qco
202	REAL p0 ! pression de reference
203	REAL p00d ! globalaverage(kpd)
204	REAL qmean_n2,qmean_n2_vert ! mass mean mixing ratio vap n2
205	REAL pqn2(ip1jmp1) ! average n2 mass index : ps*q_n2
206	REAL oldps(ip1jmp1) ! saving old pressure ps to calculate qch4
207
208	c+jld variables test conservation energie
209	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
210	C Tendance de la temp. potentiel d (theta)/ d t due a la
211	C tansformation d'energie cinetique en energie thermique
212	C cree par la dissipation
213	REAL dtetaecdt(ip1jmp1,llm)
214	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
215	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
216	! REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
217	CHARACTER*15 ztit
218	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
219	!IM SAVE ip_ebil_dyn
220	!IM DATA ip_ebil_dyn/0/
221	c-jld
222
223	! integer :: itau_w ! for write_paramLMDZ_dyn.h
224
225	! character*80 dynhist_file, dynhistave_file
226	character(len=*),parameter :: modname="leapfrog"
227	character*80 abort_message
228
229	logical dissip_conservative
230	save dissip_conservative
231	data dissip_conservative/.true./
232
233	INTEGER testita
234	PARAMETER (testita = 9)
235
236	logical , parameter :: flag_verif = .false.
237
238	! for CP(T)
239	real :: dtec
240	real :: ztetaec(ip1jmp1,llm)
241
242	! TEMP : diagnostic mass
243	real :: n2mass(iip1,jjp1)
244	real :: n2ice_ij(iip1,jjp1)
245	integer,save :: igcm_n2=0 ! index of CO2 tracer (if any)
246	integer :: i,j,ig
247	integer, parameter :: ngrid = 2+(jjm-1)*iim
248	! local mass
249	real :: mq(ip1jmp1,llm)
250
251	if (nday>=0) then
252	itaufin = nday*day_step
253	else
254	! to run a given (-nday) number of dynamical steps
255	itaufin = -nday
256	endif
257	if (less1day) then
258	c MODIF VENUS: to run less than one day:
259	itaufin = int(fractday*day_step)
260	endif
261	if (ndynstep.gt.0) then
262	! running a given number of dynamical steps
263	itaufin=ndynstep
264	endif
265	itaufinp1 = itaufin +1
266
267
268
269
270	c INITIALISATIONS
271	dudis(:,:) =0.
272	dvdis(:,:) =0.
273	dtetadis(:,:)=0.
274	dutop(:,:) =0.
275	c dvtop(:,:) =0.
276	c dtetatop(:,:)=0.
277	c dqtop(:,:,:) =0.
278	c dptop(:) =0.
279	dufi(:,:) =0.
280	dvfi(:,:) =0.
281	dtetafi(:,:)=0.
282	dqfi(:,:,:) =0.
283	dpfi(:) =0.
284
285	itau = 0
286	physics=.true.
287	if (iflag_phys==0.or.iflag_phys==2) physics=.false.
288
289	! ED18 nogcm
290	! firstcall_globaverage2d = .true.
291
292	c iday = day_ini+itau/day_step
293	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
294	c IF(time.GT.1.) THEN
295	c time = time-1.
296	c iday = iday+1
297	c ENDIF
298
299
300	c-----------------------------------------------------------------------
301	c On initialise la pression et la fonction d'Exner :
302	c --------------------------------------------------
303
304	dq(:,:,:)=0.
305	CALL pression ( ip1jmp1, ap, bp, ps, p )
306	if (pressure_exner) then
307	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
308	else
309	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
310	endif
311
312	c------------------
313	c TEST PK MONOTONE
314	c------------------
315	write(,) "Test PK"
316	do ij=1,ip1jmp1
317	do l=2,llm
318	! PRINT*,'pk(ij,l) = ',pk(ij,l)
319	! PRINT*,'pk(ij,l-1) = ',pk(ij,l-1)
320	if(pk(ij,l).gt.pk(ij,l-1)) then
321	c write(,) ij,l,pk(ij,l)
322	abort_message = 'PK non strictement decroissante'
323	call abort_gcm(modname,abort_message,1)
324	c write(,) "ATTENTION, Test PK deconnecte..."
325	endif
326	enddo
327	enddo
328	write(,) "Fin Test PK"
329	c stop
330	c------------------
331	c Preparing mixing of pressure and tracers in nogcm
332	c kpd=exp(-z/H) Needed to define a reference pressure :
333	c P0=pmean/globalaverage(kpd)
334	c thus P(i) = P0exp(-z(i)/H) = P0kpd(i)
335	c It is checked that globalaverage2d(Pi)=pmean
336	DO ij=1,ip1jmp1
337	kpd(ij) = exp(-phis(ij)/(r*200.))
338	ENDDO
339	p00d=globaverage2d(kpd) ! mean pres at ref level
340	tau_ps = 1. ! constante de rappel for pressure (s)
341	tau_n2 = 1.E5 !E5 ! constante de rappel for mix ratio qn2 (s)
342	tau_teta = 1.E7 !constante de rappel for potentiel temperature
343
344	PRINT*,'igcm_n2 = ',igcm_n2
345	do iq=1,nqtot
346	if (tname(iq)=="n2") then
347	igcm_n2=iq
348	exit
349	endif
350	enddo
351
352	c-----------------------------------------------------------------------
353	c Debut de l'integration temporelle:
354	c ----------------------------------
355
356	1 CONTINUE ! Matsuno Forward step begins here
357
358	c date: (NB: date remains unchanged for Backward step)
359	c -----
360
361	jD_cur = jD_ref + day_ini - day_ref + &
362	& (itau+1)/day_step
363	jH_cur = jH_ref + start_time + &
364	& mod(itau+1,day_step)/float(day_step)
365	jD_cur = jD_cur + int(jH_cur)
366	jH_cur = jH_cur - int(jH_cur)
367
368	c
369
370	! Save fields obtained at previous time step as '...m1'
371	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
372	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
373	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
374	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
375	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
376
377	forward = .TRUE.
378	leapf = .FALSE.
379	dt = dtvr
380
381	2 CONTINUE ! Matsuno backward or leapfrog step begins here
382
383	c-----------------------------------------------------------------------
384
385	c date: (NB: only leapfrog step requires recomputing date)
386	c -----
387
388	IF (leapf) THEN
389	jD_cur = jD_ref + day_ini - day_ref +
390	& (itau+1)/day_step
391	jH_cur = jH_ref + start_time +
392	& mod(itau+1,day_step)/float(day_step)
393	jD_cur = jD_cur + int(jH_cur)
394	jH_cur = jH_cur - int(jH_cur)
395	ENDIF
396
397
398	c gestion des appels de la physique et des dissipations:
399	c ------------------------------------------------------
400	c
401	c ... P.Le Van ( 6/02/95 ) ....
402
403	! ED18: suppression des mentions de la variable apdiss dans le cas
404	! 'nogcm'
405
406	apphys = .FALSE.
407	statcl = .FALSE.
408	conser = .FALSE.
409
410
411	IF( purmats ) THEN
412	! Purely Matsuno time stepping
413	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
414
415
416	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
417	s .and. physics ) apphys = .TRUE.
418	ELSE
419	! Leapfrog/Matsuno time stepping
420	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
421
422
423	IF( MOD(itau+1,iphysiq).EQ.0.AND.physics) apphys=.TRUE.
424	END IF
425
426	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
427	! supress dissipation step
428
429
430
431	c-----------------------------------------------------------------------
432	c calcul des tendances dynamiques:
433	c --------------------------------
434
435	! ED18: suppression de l'onglet pour le cas nogcm
436
437
438	dv(:,:) = 0.
439	du(:,:) = 0.
440	dteta(:,:) = 0.
441	dq(:,:,:) = 0.
442
443
444
445	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
446	c
447	c-----------------------------------------------------------------------
448	c calcul des tendances physiques:
449	c -------------------------------
450	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
451	c
452	IF( purmats ) THEN
453	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
454	ELSE
455	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
456	ENDIF
457	c
458	c
459	IF( apphys ) THEN
460	c
461	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
462	c
463
464	CALL pression ( ip1jmp1, ap, bp, ps, p )
465	if (pressure_exner) then
466	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
467	else
468	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
469	endif
470
471	! Compute geopotential (physics might need it)
472	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
473
474	jD_cur = jD_ref + day_ini - day_ref + &
475	& (itau+1)/day_step
476
477	! AS: we make jD_cur to be pday
478	jD_cur = int(day_ini + itau/day_step)
479
480	! print*,'itau =',itau
481	! print*,'day_step =',day_step
482	! print*,'itau/day_step =',itau/day_step
483
484
485	jH_cur = jH_ref + start_time + &
486	& mod(itau+1,day_step)/float(day_step)
487	jH_cur = jH_ref + start_time + &
488	& mod(itau,day_step)/float(day_step)
489	jD_cur = jD_cur + int(jH_cur)
490	jH_cur = jH_cur - int(jH_cur)
491
492
493
494	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
495	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
496	! write(lunout,*)'current date = ',an, mois, jour, secondes
497
498	c rajout debug
499	c lafin = .true.
500
501
502	c Interface avec les routines de phylmd (phymars ... )
503	c -----------------------------------------------------
504
505	c+jld
506
507	c Diagnostique de conservation de l'Energie : initialisation
508	IF (ip_ebil_dyn.ge.1 ) THEN
509	ztit='bil dyn'
510	! Ehouarn: be careful, diagedyn is Earth-specific!
511	IF (planet_type.eq."earth") THEN
512	CALL diagedyn(ztit,2,1,1,dtphys
513	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
514	ENDIF
515	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
516	c-jld
517	#ifdef CPP_IOIPSL
518	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
519	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
520	c IF (first) THEN
521	c first=.false.
522	c#include "ini_paramLMDZ_dyn.h"
523	c ENDIF
524	c
525	c#include "write_paramLMDZ_dyn.h"
526	c
527	#endif
528	! #endif of #ifdef CPP_IOIPSL
529
530	c call WriteField('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
531	! print*,'---------LEAPFROG---------------'
532	! print*,''
533	! print*,'> AVANT CALFIS'
534	! print*,''
535	! print*,'teta(3049,:) = ',teta(3049,:)
536	! print*,''
537	! print*,'dteta(3049,:) = ',dteta(3049,:)
538	! print*,''
539	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
540	! print*,''
541
542
543	CALL calfis( lafin , jD_cur, jH_cur,
544	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
545	$ du,dv,dteta,dq,
546	$ flxw,
547	$ dufi,dvfi,dtetafi,dqfi,dpfi )
548
549	c call WriteField('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
550	c call WriteField('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
551	c call WriteField('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
552	! print*,'> APRES CALFIS (AVANT ADDFI)'
553	! print*,''
554	! print*,'teta(3049,:) = ',teta(3049,:)
555	! print*,''
556	! print*,'dteta(3049,:) = ',dteta(3049,:)
557	! print*,''
558	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
559	! print*,''
560
561
562	c ajout des tendances physiques
563	c ------------------------------
564
565	CALL addfi( dtphys, leapf, forward ,
566	$ ucov, vcov, teta , q ,ps ,
567	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
568
569	! print*,'> APRES ADDFI'
570	! print*,''
571	! print*,'teta(3049,:) = ',teta(3049,:)
572	! print*,''
573	! print*,'dteta(3049,:) = ',dteta(3049,:)
574	! print*,''
575	! print*,'dtetafi(3049,:) = ',dtetafi(3049,:)
576	! print*,''
577
578	CALL pression (ip1jmp1,ap,bp,ps,p)
579	CALL massdair(p,masse)
580
581	! Mass of tracers in each mess (kg) before Horizontal mixing of pressure
582	c -------------------------------------------------------------------
583
584	DO l=1, llm
585	DO ij=1,ip1jmp1
586	mq(ij,l) = masse(ij,l)*q(ij,l,igcm_n2)
587	ENDDO
588	ENDDO
589
590	c Horizontal mixing of pressure
591	c ------------------------------
592	! Rappel newtonien vers psmean
593	psmean= globaverage2d(ps) ! mean pressure
594	p0=psmean/p00d
595	DO ij=1,ip1jmp1
596	oldps(ij)=ps(ij)
597	ps(ij)=ps(ij) +(p0*kpd(ij)
598	& -ps(ij))*(1-exp(-dtphys/tau_ps))
599	ENDDO
600
601	write(72,*) 'psmean = ',psmean ! mean pressure redistributed
602
603	! security check: test pressure negative
604	DO ij=1,ip1jmp1
605	IF (ps(ij).le.0.) then
606	!write(,) 'Negative pressure :'
607	!write(,) 'ps = ',ps(ij),' ig = ',ij
608	!write(,) 'pmean = ',p0*kpd(ij)
609	!write(,) 'tau_ps = ',tau_ps
610	!STOP
611	ps(ij)=0.0000001*kpd(ij)/p00d
612	ENDIF
613	ENDDO
614	!***********************
615	! Correction on teta due to surface pressure changes
616	DO l = 1,llm
617	DO ij = 1,ip1jmp1
618	teta(ij,l)= teta(ij,l)(ps(ij)/oldps(ij))*kappa
619	ENDDO
620	ENDDO
621	!***********************
622
623
624
625	! ! update pressure and update p and pk
626	! DO ij=1,ip1jmp1
627	! ps(ij) = ps(ij) + dpfi(ij)*dtphys
628	! ENDDO
629	CALL pression (ip1jmp1,ap,bp,ps,p)
630	CALL massdair(p,masse)
631	if (pressure_exner) then
632	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
633	else
634	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
635	endif
636
637	! Update tracers after Horizontal mixing of pressure to ! conserve tracer mass
638	c -------------------------------------------------------------------
639	DO l=1, llm
640	DO ij=1,ip1jmp1
641	q(ij,l,igcm_n2) = mq(ij,l)/ masse(ij,l)
642	ENDDO
643	ENDDO
644
645
646
647	c Horizontal mixing of pressure
648	c ------------------------------
649	! Rappel newtonien vers psmean
650	psmean= globaverage2d(ps) ! mean pressure
651	! ! increment q_n2 with physical tendancy
652	! IF (igcm_n2.ne.0) then
653	! DO l=1, llm
654	! DO ij=1,ip1jmp1
655	! q(ij,l,igcm_n2)=q(ij,l,igcm_n2)+
656	! & dqfi(ij,l,igcm_n2)*dtphys
657	! ENDDO
658	! ENDDO
659	! ENDIF
660
661	c Mixing CO2 vertically
662	c --------------------------
663	if (igcm_n2.ne.0) then
664	DO ij=1,ip1jmp1
665	qmean_n2_vert=0.
666	DO l=1, llm
667	qmean_n2_vert= qmean_n2_vert
668	& + q(ij,l,igcm_n2)*( p(ij,l) - p(ij,l+1))
669	END DO
670	qmean_n2_vert= qmean_n2_vert/ps(ij)
671	DO l=1, llm
672	q(ij,l,igcm_n2)= qmean_n2_vert
673	END DO
674	END DO
675	end if
676
677
678
679	c Horizontal mixing of pressure
680	c ------------------------------
681	! Rappel newtonien vers psmean
682	psmean= globaverage2d(ps) ! mean pressure
683
684	c Horizontal mixing tracers, and temperature (replace dynamics in nogcm)
685	c ---------------------------------------------------------------------------------
686
687	! Simulate redistribution by dynamics for qn2
688	if (igcm_n2.ne.0) then
689
690	DO ij=1,ip1jmp1
691	pqn2(ij)= ps(ij) * q(ij,1,igcm_n2)
692	ENDDO
693	pqn2mean=globaverage2d(pqn2)
694
695	! Rappel newtonien vers qn2_mean
696	qmean_n2= pqn2mean / psmean
697
698	DO ij=1,ip1jmp1
699	q(ij,1,igcm_n2)=q(ij,1,igcm_n2)+
700	& (qmean_n2-q(ij,1,igcm_n2))*
701	& (1.-exp(-dtphys/tau_n2))
702	ENDDO
703
704	DO l=2, llm
705	DO ij=1,ip1jmp1
706	q(ij,l,igcm_n2)=q(ij,1,igcm_n2)
707	END DO
708	END DO
709
710	! TEMPORAIRE (ED)
711	! PRINT*,'psmean = ',psmean
712	! PRINT*,'qmean_n2 = ',qmean_n2
713	! PRINT*,'pqn2mean = ',pqn2mean
714	! PRINT*,'q(50,1,igcm_n2) = ',q(50,1,igcm_n2)
715	! PRINT*,'q(50,2,igcm_n2) = ',q(50,2,igcm_n2)
716	! PRINT*,'q(50,3,igcm_n2) = ',q(50,3,igcm_n2)
717
718	endif ! igcm_n2.ne.0
719
720
721	! ********************************************************
722
723	c Horizontal mixing of pressure
724	c ------------------------------
725	! Rappel newtonien vers psmean
726	psmean= globaverage2d(ps) ! mean pressure
727
728
729	c Mixing Temperature horizontally
730	c -------------------------------
731	! initialize variables that will be averaged
732	DO l=1,llm
733	DO ij=1,ip1jmp1
734	dp(ij,l) = p(ij,l) - p(ij,l+1)
735	tetadp(ij,l) = teta(ij,l)*dp(ij,l)
736	ENDDO
737	ENDDO
738
739	DO l=1,llm
740	tetadpmean = globaverage2d(tetadp(:,l))
741	dpmean = globaverage2d(dp(:,l))
742	tetamean = tetadpmean / dpmean
743	DO ij=1,ip1jmp1
744	teta(ij,l) = teta(ij,l) + (tetamean - teta(ij,l)) *
745	& (1 - exp(-dtphys/tau_teta))
746	ENDDO
747	ENDDO
748
749
750	ENDIF ! of IF( apphys )
751
752
753	c ********************************************************************
754	c ********************************************************************
755	c .... fin de l'integration dynamique et physique pour le pas itau ..
756	c ********************************************************************
757	c ********************************************************************
758
759
760	IF ( .NOT.purmats ) THEN
761	c ........................................................
762	c .............. schema matsuno + leapfrog ..............
763	c ........................................................
764
765	IF(forward. OR. leapf) THEN
766	itau= itau + 1
767	c iday= day_ini+itau/day_step
768	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
769	c IF(time.GT.1.) THEN
770	c time = time-1.
771	c iday = iday+1
772	c ENDIF
773	ENDIF
774
775
776	IF( itau. EQ. itaufinp1 ) then
777	if (flag_verif) then
778	write(79,*) 'ucov',ucov
779	write(80,*) 'vcov',vcov
780	write(81,*) 'teta',teta
781	write(82,*) 'ps',ps
782	write(83,*) 'q',q
783	WRITE(85,*) 'q1 = ',q(:,:,1)
784	WRITE(86,*) 'q3 = ',q(:,:,3)
785	endif
786
787	abort_message = 'Simulation finished'
788
789	call abort_gcm(modname,abort_message,0)
790	ENDIF
791	c-----------------------------------------------------------------------
792	c ecriture du fichier histoire moyenne:
793	c -------------------------------------
794
795	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
796	IF(itau.EQ.itaufin) THEN
797	iav=1
798	ELSE
799	iav=0
800	ENDIF
801
802	! ! Ehouarn: re-compute geopotential for outputs
803	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
804
805	IF (ok_dynzon) THEN
806	#ifdef CPP_IOIPSL
807	c les traceurs ne sont pas sortis, trop lourd.
808	c Peut changer eventuellement si besoin.
809	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
810	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
811	& du,dudis,dutop,dufi)
812	#endif
813	END IF
814	IF (ok_dyn_ave) THEN
815	#ifdef CPP_IOIPSL
816	CALL writedynav(itau,vcov,
817	& ucov,teta,pk,phi,q,masse,ps,phis)
818	#endif
819	ENDIF
820
821	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
822
823	if (ok_iso_verif) then
824	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
825	endif !if (ok_iso_verif) then
826
827	c-----------------------------------------------------------------------
828	c ecriture de la bande histoire:
829	c ------------------------------
830
831	IF( MOD(itau,iecri).EQ.0) THEN
832	! Ehouarn: output only during LF or Backward Matsuno
833	if (leapf.or.(.not.leapf.and.(.not.forward))) then
834	! ADAPTATION GCM POUR CP(T)
835	call tpot2t(ijp1llm,teta,temp,pk)
836	tsurpk = cpp*temp/pk
837	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
838	unat=0.
839	do l=1,llm
840	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
841	vnat(:,l)=vcov(:,l)/cv(:)
842	enddo
843	#ifdef CPP_IOIPSL
844	if (ok_dyn_ins) then
845	! write(lunout,*) "leapfrog: call writehist, itau=",itau
846	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
847	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
848	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
849	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
850	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
851	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
852	endif ! of if (ok_dyn_ins)
853	#endif
854	! For some Grads outputs of fields
855	if (output_grads_dyn) then
856	#include "write_grads_dyn.h"
857	endif
858	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
859	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
860
861	IF(itau.EQ.itaufin) THEN
862
863	CALL dynredem1("restart.nc",start_time,
864	& vcov,ucov,teta,q,masse,ps)
865	CLOSE(99)
866	!!! Ehouarn: Why not stop here and now?
867	ENDIF ! of IF (itau.EQ.itaufin)
868
869	c-----------------------------------------------------------------------
870	c gestion de l'integration temporelle:
871	c ------------------------------------
872
873	IF( MOD(itau,iperiod).EQ.0 ) THEN
874	GO TO 1
875	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
876
877	IF( forward ) THEN
878	c fin du pas forward et debut du pas backward
879
880	forward = .FALSE.
881	leapf = .FALSE.
882	GO TO 2
883
884	ELSE
885	c fin du pas backward et debut du premier pas leapfrog
886
887	leapf = .TRUE.
888	dt = 2.*dtvr
889	GO TO 2
890	END IF ! of IF (forward)
891	ELSE
892
893	c ...... pas leapfrog .....
894
895	leapf = .TRUE.
896	dt = 2.*dtvr
897	GO TO 2
898	END IF ! of IF (MOD(itau,iperiod).EQ.0)
899	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
900
901	ELSE ! of IF (.not.purmats)
902
903	if (ok_iso_verif) then
904	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
905	endif !if (ok_iso_verif) then
906
907	c ........................................................
908	c .............. schema matsuno ...............
909	c ........................................................
910	IF( forward ) THEN
911
912	itau = itau + 1
913	c iday = day_ini+itau/day_step
914	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
915	c
916	c IF(time.GT.1.) THEN
917	c time = time-1.
918	c iday = iday+1
919	c ENDIF
920
921	forward = .FALSE.
922	IF( itau. EQ. itaufinp1 ) then
923	abort_message = 'Simulation finished'
924	call abort_gcm(modname,abort_message,0)
925	ENDIF
926	GO TO 2
927
928	ELSE ! of IF(forward) i.e. backward step
929
930	if (ok_iso_verif) then
931	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
932	endif !if (ok_iso_verif) then
933
934	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
935	IF(itau.EQ.itaufin) THEN
936	iav=1
937	ELSE
938	iav=0
939	ENDIF
940
941	! ! Ehouarn: re-compute geopotential for outputs
942	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
943
944	IF (ok_dynzon) THEN
945	#ifdef CPP_IOIPSL
946	c les traceurs ne sont pas sortis, trop lourd.
947	c Peut changer eventuellement si besoin.
948	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
949	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
950	& du,dudis,dutop,dufi)
951	#endif
952	ENDIF
953	IF (ok_dyn_ave) THEN
954	#ifdef CPP_IOIPSL
955	CALL writedynav(itau,vcov,
956	& ucov,teta,pk,phi,q,masse,ps,phis)
957	#endif
958	ENDIF
959
960	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
961
962	IF(MOD(itau,iecri ).EQ.0) THEN
963	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
964	! ADAPTATION GCM POUR CP(T)
965	call tpot2t(ijp1llm,teta,temp,pk)
966	tsurpk = cpp*temp/pk
967	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
968	unat=0.
969	do l=1,llm
970	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
971	vnat(:,l)=vcov(:,l)/cv(:)
972	enddo
973	#ifdef CPP_IOIPSL
974	if (ok_dyn_ins) then
975	! write(lunout,*) "leapfrog: call writehist (b)",
976	! & itau,iecri
977	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
978	endif ! of if (ok_dyn_ins)
979	#endif
980	! For some Grads outputs
981	if (output_grads_dyn) then
982	#include "write_grads_dyn.h"
983	endif
984
985	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
986
987	IF(itau.EQ.itaufin) THEN
988	CALL dynredem1("restart.nc",start_time,
989	& vcov,ucov,teta,q,masse,ps)
990	ENDIF ! of IF(itau.EQ.itaufin)
991
992	forward = .TRUE.
993	GO TO 1
994
995	ENDIF ! of IF (forward)
996
997	END IF ! of IF(.not.purmats)
998
999	STOP
1000	END
1001
1002	! ************************************************************
1003	! globalaverage VERSION PLuto
1004	! ************************************************************
1005	FUNCTION globaverage2d(var)
1006	! ************************************************************
1007	IMPLICIT NONE
1008	#include "dimensions.h"
1009	#include "paramet.h"
1010	#include "comgeom2.h"
1011	!ip1jmp1 called in paramet.h = iip1 x jjp1
1012	REAL var(iip1,jjp1), globaverage2d
1013	integer i,j
1014	REAL airetot
1015	save airetot
1016	logical firstcall
1017	data firstcall/.true./
1018	save firstcall
1019
1020	if (firstcall) then
1021	firstcall=.false.
1022	airetot =0.
1023	DO j=2,jjp1-1 ! lat
1024	DO i = 1, iim ! lon
1025	airetot = airetot + aire(i,j)
1026	END DO
1027	END DO
1028	DO i=1,iim
1029	airetot=airetot + aire(i,1)+aire(i,jjp1)
1030	ENDDO
1031	end if
1032
1033	globaverage2d = 0.
1034	DO j=2,jjp1-1
1035	DO i = 1, iim
1036	globaverage2d = globaverage2d + var(i,j)*aire(i,j)
1037	END DO
1038	END DO
1039
1040	DO i=1,iim
1041	globaverage2d = globaverage2d + var(i,1)*aire(i,1)
1042	globaverage2d = globaverage2d + var(i,jjp1)*aire(i,jjp1)
1043	ENDDO
1044
1045	globaverage2d = globaverage2d/airetot
1046	return
1047	end
1048

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

Download in other formats:

Original Format