Context Navigation

iniacademic.F @ 3338

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

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

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

File size: 8.4 KB

Line
1	!
2	! $Id: iniacademic.F 1446 2010-10-22 09:27:25Z acozic $
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, ONLY: day_step,planet_type
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	c
98	! initialize planet radius, rotation rate,...
99	call conf_planete
100
101	time_0=0.
102	day_ref=1
103	annee_ref=0
104
105	im = iim
106	jm = jjm
107	day_ini = 1
108	dtvr = daysec/REAL(day_step)
109	zdtvr=dtvr
110	etot0 = 0.
111	ptot0 = 0.
112	ztot0 = 0.
113	stot0 = 0.
114	ang0 = 0.
115
116	if (llm.eq.1) then
117	! specific initializations for the shallow water case
118	kappa=1
119	endif
120
121	CALL iniconst
122	CALL inigeom
123	CALL inifilr
124
125	if (llm.eq.1) then
126	! initialize fields for the shallow water case, if required
127	if (.not.read_start) then
128	phis(:)=0.
129	q(:,:,:)=0
130	CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
131	endif
132	endif
133
134	if (iflag_phys.eq.2) then
135	! initializations for the academic case
136
137	! if (planet_type=="earth") then
138
139	! 1. local parameters
140	! by convention, winter is in the southern hemisphere
141	! Geostrophic wind or no wind?
142	ok_geost=.TRUE.
143	CALL getin('ok_geost',ok_geost)
144	! Constants for Newtonian relaxation and friction
145	k_f=1. !friction
146	CALL getin('k_j',k_f)
147	k_f=1./(daysec*k_f)
148	k_c_s=4. !cooling surface
149	CALL getin('k_c_s',k_c_s)
150	k_c_s=1./(daysec*k_c_s)
151	k_c_a=40. !cooling free atm
152	CALL getin('k_c_a',k_c_a)
153	k_c_a=1./(daysec*k_c_a)
154	! Constants for Teta equilibrium profile
155	teta0=315. ! mean Teta (S.H. 315K)
156	CALL getin('teta0',teta0)
157	ttp=200. ! Tropopause temperature (S.H. 200K)
158	CALL getin('ttp',ttp)
159	eps=0. ! Deviation to N-S symmetry(~0-20K)
160	CALL getin('eps',eps)
161	delt_y=60. ! Merid Temp. Gradient (S.H. 60K)
162	CALL getin('delt_y',delt_y)
163	delt_z=10. ! Vertical Gradient (S.H. 10K)
164	CALL getin('delt_z',delt_z)
165	! Polar vortex
166	ok_pv=.false.
167	CALL getin('ok_pv',ok_pv)
168	phi_pv=-50. ! Latitude of edge of vortex
169	CALL getin('phi_pv',phi_pv)
170	phi_pv=phi_pv*pi/180.
171	dphi_pv=5. ! Width of the edge
172	CALL getin('dphi_pv',dphi_pv)
173	dphi_pv=dphi_pv*pi/180.
174	gam_pv=4. ! -dT/dz vortex (in K/km)
175	CALL getin('gam_pv',gam_pv)
176
177	! 2. Initialize fields towards which to relax
178	! Friction
179	knewt_g=k_c_a
180	DO l=1,llm
181	zsig=presnivs(l)/preff
182	knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
183	kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
184	ENDDO
185	DO j=1,jjp1
186	clat4((j-1)iip1+1:jiip1)=cos(rlatu(j))**4
187	ENDDO
188
189	! Potential temperature
190	DO l=1,llm
191	zsig=presnivs(l)/preff
192	tetastrat=ttpzsig*(-kappa)
193	tetapv=tetastrat
194	IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
195	tetapv=tetastrat(zsig10.)*(kappacpp*gam_pv/1000./g)
196	ENDIF
197	DO j=1,jjp1
198	! Troposphere
199	ddsin=sin(rlatu(j))
200	tetajl(j,l)=teta0-delt_yddsinddsin+eps*ddsin
201	& -delt_z(1.-ddsinddsin)*log(zsig)
202	! Profil stratospherique isotherme (+vortex)
203	w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
204	tetastrat=tetastrat(1.-w_pv)+tetapvw_pv
205	tetajl(j,l)=MAX(tetajl(j,l),tetastrat)
206	ENDDO
207	ENDDO ! of DO l=1,llm
208
209	! CALL writefield('theta_eq',tetajl)
210
211	do l=1,llm
212	do j=1,jjp1
213	do i=1,iip1
214	ij=(j-1)*iip1+i
215	tetarappel(ij,l)=tetajl(j,l)
216	enddo
217	enddo
218	enddo
219
220
221	! else
222	! write(lunout,*)"iniacademic: planet types other than earth",
223	! & " not implemented (yet)."
224	! stop
225	! endif ! of if (planet_type=="earth")
226
227	! 3. Initialize fields (if necessary)
228	IF (.NOT. read_start) THEN
229	! surface pressure
230	ps(:)=preff
231	! ground geopotential
232	phis(:)=0.
233
234	CALL pression ( ip1jmp1, ap, bp, ps, p )
235	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
236	CALL massdair(p,masse)
237
238	! bulk initialization of temperature
239	teta(:,:)=tetarappel(:,:)
240
241	! geopotential
242	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
243
244	! winds
245	if (ok_geost) then
246	call ugeostr(phi,ucov)
247	else
248	ucov(:,:)=0.
249	endif
250	vcov(:,:)=0.
251
252	! bulk initialization of tracers
253	if (planet_type=="earth") then
254	! Earth: first two tracers will be water
255	do i=1,nqtot
256	if (i.eq.1) q(:,:,i)=1.e-10
257	if (i.eq.2) q(:,:,i)=1.e-15
258	if (i.gt.2) q(:,:,i)=0.
259	enddo
260	else
261	q(:,:,:)=0
262	endif ! of if (planet_type=="earth")
263
264	! add random perturbation to temperature
265	idum = -1
266	zz = ran1(idum)
267	idum = 0
268	do l=1,llm
269	do ij=iip2,ip1jm
270	teta(ij,l)=teta(ij,l)(1.+0.005ran1(idum))
271	enddo
272	enddo
273
274	! maintain periodicity in longitude
275	do l=1,llm
276	do ij=1,ip1jmp1,iip1
277	teta(ij+iim,l)=teta(ij,l)
278	enddo
279	enddo
280
281	ENDIF ! of IF (.NOT. read_start)
282	endif ! of if (iflag_phys.eq.2)
283
284	END
285	c-----------------------------------------------------------------------

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

Download in other formats:

Original Format