Context Navigation

iniacademic.F @ 1437

Last change on this file since 1437 was 1437, checked in by Ehouarn Millour, 14 years ago

Implementation of FC's improved Newtonian configuration: more realistic (and tunable) friction & fields towards which to relax and planet parameters (gravity, radius, rotation rate...) can be read in from a planet.def file (this occurs only when running in "newtonian mode"; planet parameters are otherwise read in from the start file).

Checked that these implementation don't change results when running standard Earth bench tests on Vargas.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 8.5 KB

Line
1	!
2	! $Id: iniacademic.F 1437 2010-09-30 08:29:10Z emillour $
3	!
4	c
5	c
6	SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
7
8	USE filtreg_mod
9	USE infotrac, ONLY : nqtot
10	USE control_mod
11	#ifdef CPP_IOIPSL
12	USE IOIPSL
13	#else
14	! if not using IOIPSL, we still need to use (a local version of) getin
15	USE ioipsl_getincom
16	#endif
17	USE Write_Field
18
19
20	c%W% %G%
21	c=======================================================================
22	c
23	c Author: Frederic Hourdin original: 15/01/93
24	c -------
25	c
26	c Subject:
27	c ------
28	c
29	c Method:
30	c --------
31	c
32	c Interface:
33	c ----------
34	c
35	c Input:
36	c ------
37	c
38	c Output:
39	c -------
40	c
41	c=======================================================================
42	IMPLICIT NONE
43	c-----------------------------------------------------------------------
44	c Declararations:
45	c ---------------
46
47	#include "dimensions.h"
48	#include "paramet.h"
49	#include "comvert.h"
50	#include "comconst.h"
51	#include "comgeom.h"
52	#include "academic.h"
53	#include "ener.h"
54	#include "temps.h"
55	#include "iniprint.h"
56	#include "logic.h"
57
58	c Arguments:
59	c ----------
60
61	real time_0
62
63	c variables dynamiques
64	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
65	REAL teta(ip1jmp1,llm) ! temperature potentielle
66	REAL q(ip1jmp1,llm,nqtot) ! champs advectes
67	REAL ps(ip1jmp1) ! pression au sol
68	REAL masse(ip1jmp1,llm) ! masse d'air
69	REAL phis(ip1jmp1) ! geopotentiel au sol
70
71	c Local:
72	c ------
73
74	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
75	REAL pks(ip1jmp1) ! exner au sol
76	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
77	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
78	REAL phi(ip1jmp1,llm) ! geopotentiel
79	REAL ddsin,tetastrat,zsig,tetapv,w_pv ! variables auxiliaires
80	real tetajl(jjp1,llm)
81	INTEGER i,j,l,lsup,ij
82
83	REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
84	REAL k_f,k_c_a,k_c_s ! Constantes de rappel
85	LOGICAL ok_geost ! Initialisation vent geost. ou nul
86	LOGICAL ok_pv ! Polar Vortex
87	REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex
88
89	real zz,ran1
90	integer idum
91
92	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
93
94	c-----------------------------------------------------------------------
95	! 1. Initializations for Earth-like case
96	! --------------------------------------
97	if (planet_type=="earth") then
98	c
99	! initialize planet radius, rotation rate,...
100	call conf_planete
101
102	time_0=0.
103	day_ref=1
104	annee_ref=0
105
106	im = iim
107	jm = jjm
108	day_ini = 1
109	dtvr = daysec/REAL(day_step)
110	zdtvr=dtvr
111	etot0 = 0.
112	ptot0 = 0.
113	ztot0 = 0.
114	stot0 = 0.
115	ang0 = 0.
116
117	if (llm.eq.1) then
118	! specific initializations for the shallow water case
119	kappa=1
120	endif
121
122	CALL iniconst
123	CALL inigeom
124	CALL inifilr
125
126	if (llm.eq.1) then
127	! initialize fields for the shallow water case, if required
128	if (.not.read_start) then
129	phis(:)=0.
130	q(:,:,1)=1.e-10
131	q(:,:,2)=1.e-15
132	q(:,:,3:nqtot)=0.
133	CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
134	endif
135	endif
136
137	if (iflag_phys.eq.2) then
138	! initializations for the academic case
139
140	! 1. local parameters
141	! by convention, winter is in the southern hemisphere
142	! Geostrophic wind or no wind?
143	ok_geost=.TRUE.
144	CALL getin('ok_geost',ok_geost)
145	! Constants for Newtonian relaxation and friction
146	k_f=1. !friction
147	CALL getin('k_j',k_f)
148	k_f=1./(daysec*k_f)
149	k_c_s=4. !cooling surface
150	CALL getin('k_c_s',k_c_s)
151	k_c_s=1./(daysec*k_c_s)
152	k_c_a=40. !cooling free atm
153	CALL getin('k_c_a',k_c_a)
154	k_c_a=1./(daysec*k_c_a)
155	! Constants for Teta equilibrium profile
156	teta0=315. ! mean Teta (S.H. 315K)
157	CALL getin('teta0',teta0)
158	ttp=200. ! Tropopause temperature (S.H. 200K)
159	CALL getin('ttp',ttp)
160	eps=0. ! Deviation to N-S symmetry(~0-20K)
161	CALL getin('eps',eps)
162	delt_y=60. ! Merid Temp. Gradient (S.H. 60K)
163	CALL getin('delt_y',delt_y)
164	delt_z=10. ! Vertical Gradient (S.H. 10K)
165	CALL getin('delt_z',delt_z)
166	! Polar vortex
167	ok_pv=.false.
168	CALL getin('ok_pv',ok_pv)
169	phi_pv=-50. ! Latitude of edge of vortex
170	CALL getin('phi_pv',phi_pv)
171	phi_pv=phi_pv*pi/180.
172	dphi_pv=5. ! Width of the edge
173	CALL getin('dphi_pv',dphi_pv)
174	dphi_pv=dphi_pv*pi/180.
175	gam_pv=4. ! -dT/dz vortex (in K/km)
176	CALL getin('gam_pv',gam_pv)
177
178	! 2. Initialize fields towards which to relax
179	! Friction
180	knewt_g=k_c_a
181	DO l=1,llm
182	zsig=presnivs(l)/preff
183	knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
184	kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
185	ENDDO
186	DO j=1,jjp1
187	clat4((j-1)iip1+1:jiip1)=cos(rlatu(j))**4
188	ENDDO
189
190	! Potential temperature
191	DO l=1,llm
192	zsig=presnivs(l)/preff
193	tetastrat=ttpzsig*(-kappa)
194	tetapv=tetastrat
195	IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
196	tetapv=tetastrat(zsig10.)*(kappacpp*gam_pv/1000./g)
197	ENDIF
198	DO j=1,jjp1
199	! Troposphere
200	ddsin=sin(rlatu(j))
201	tetajl(j,l)=teta0-delt_yddsinddsin+eps*ddsin
202	& -delt_z(1.-ddsinddsin)*log(zsig)
203	! Profil stratospherique isotherme (+vortex)
204	w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
205	tetastrat=tetastrat(1.-w_pv)+tetapvw_pv
206	tetajl(j,l)=MAX(tetajl(j,l),tetastrat)
207	ENDDO
208	ENDDO ! of DO l=1,llm
209
210	! CALL writefield('theta_eq',tetajl)
211
212	do l=1,llm
213	do j=1,jjp1
214	do i=1,iip1
215	ij=(j-1)*iip1+i
216	tetarappel(ij,l)=tetajl(j,l)
217	enddo
218	enddo
219	enddo
220
221	! 3. Initialize fields (if necessary)
222	IF (.NOT. read_start) THEN
223	! surface pressure
224	ps(:)=preff
225	! ground geopotential
226	phis(:)=0.
227
228	CALL pression ( ip1jmp1, ap, bp, ps, p )
229	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
230	CALL massdair(p,masse)
231
232	! bulk initialization of temperature
233	teta(:,:)=tetarappel(:,:)
234
235	! geopotential
236	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
237
238	! winds
239	if (ok_geost) then
240	call ugeostr(phi,ucov)
241	else
242	ucov(:,:)=0.
243	endif
244	vcov(:,:)=0.
245
246	! bulk initialization of tracers
247	do i=1,nqtot
248	if (i.eq.1) q(:,:,i)=1.e-10
249	if (i.eq.2) q(:,:,i)=1.e-15
250	if (i.gt.2) q(:,:,i)=0.
251	enddo
252
253	! add random perturbation to temperature
254	idum = -1
255	zz = ran1(idum)
256	idum = 0
257	do l=1,llm
258	do ij=iip2,ip1jm
259	teta(ij,l)=teta(ij,l)(1.+0.005ran1(idum))
260	enddo
261	enddo
262
263	do l=1,llm
264	do ij=1,ip1jmp1,iip1
265	teta(ij+iim,l)=teta(ij,l)
266	enddo
267	enddo
268
269	c PRINT *,' Appel test_period avec tetarappel '
270	c CALL test_period ( ucov,vcov,tetarappel,q,p,phis )
271	c PRINT *,' Appel test_period avec teta '
272	c CALL test_period ( ucov,vcov,teta,q,p,phis )
273
274	! initialize a traceur on one column
275	! j=jjp1*3/4
276	! i=iip1/2
277	! ij=(j-1)*iip1+i
278	! q(ij,:,3)=1.
279
280	ENDIF ! of IF (.NOT. read_start)
281	endif ! of if (iflag_phys.eq.2)
282
283	else
284	write(lunout,*)"iniacademic: planet types other than earth",
285	& " not implemented (yet)."
286	stop
287	endif ! of if (planet_type=="earth")
288	return
289	END
290	c-----------------------------------------------------------------------

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

Download in other formats:

Original Format