Context Navigation

iniacademic_loc.F90 @ 2253

Last change on this file since 2253 was 2087, checked in by fhourdin, 10 years ago
Correction de bug pour gfortran qui refuse if (rerad_start==.false.) change en if (.not. read_start)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 9.5 KB

Line
1	!
2	! $Id: iniacademic.F90 1625 2012-05-09 13:14:48Z lguez $
3	!
4	SUBROUTINE iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
5
6	USE filtreg_mod, ONLY: inifilr
7	use exner_hyb_m, only: exner_hyb
8	use exner_milieu_m, only: exner_milieu
9	USE infotrac, ONLY : nqtot
10	USE control_mod, ONLY: day_step,planet_type
11	USE parallel_lmdz, ONLY: ijb_u, ije_u, ijb_v, ije_v
12	#ifdef CPP_IOIPSL
13	USE IOIPSL, ONLY: getin
14	#else
15	! if not using IOIPSL, we still need to use (a local version of) getin
16	USE ioipsl_getincom, ONLY: getin
17	#endif
18	USE Write_Field
19
20	! Author: Frederic Hourdin original: 15/01/93
21	! The forcing defined here is from Held and Suarez, 1994, Bulletin
22	! of the American Meteorological Society, 75, 1825.
23
24	IMPLICIT NONE
25
26	! Declararations:
27	! ---------------
28
29	include "dimensions.h"
30	include "paramet.h"
31	include "comvert.h"
32	include "comconst.h"
33	include "comgeom.h"
34	include "academic.h"
35	include "ener.h"
36	include "temps.h"
37	include "iniprint.h"
38	include "logic.h"
39
40	! Arguments:
41	! ----------
42
43	REAL,INTENT(OUT) :: time_0
44
45	! fields
46	REAL,INTENT(OUT) :: vcov(ijb_v:ije_v,llm) ! meridional covariant wind
47	REAL,INTENT(OUT) :: ucov(ijb_u:ije_u,llm) ! zonal covariant wind
48	REAL,INTENT(OUT) :: teta(ijb_u:ije_u,llm) ! potential temperature (K)
49	REAL,INTENT(OUT) :: q(ijb_u:ije_u,llm,nqtot) ! advected tracers (.../kg_of_air)
50	REAL,INTENT(OUT) :: ps(ijb_u:ije_u) ! surface pressure (Pa)
51	REAL,INTENT(OUT) :: masse(ijb_u:ije_u,llm) ! air mass in grid cell (kg)
52	REAL,INTENT(OUT) :: phis(ijb_u:ije_u) ! surface geopotential
53
54	! Local:
55	! ------
56
57	REAL,ALLOCATABLE :: vcov_glo(:,:),ucov_glo(:,:),teta_glo(:,:)
58	REAL,ALLOCATABLE :: q_glo(:,:),masse_glo(:,:),ps_glo(:)
59	REAL,ALLOCATABLE :: phis_glo(:)
60	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
61	REAL pks(ip1jmp1) ! exner au sol
62	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
63	REAL phi(ip1jmp1,llm) ! geopotentiel
64	REAL ddsin,zsig,tetapv,w_pv ! variables auxiliaires
65	real tetastrat ! potential temperature in the stratosphere, in K
66	real tetajl(jjp1,llm)
67	INTEGER i,j,l,lsup,ij
68
69	REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
70	REAL k_f,k_c_a,k_c_s ! Constantes de rappel
71	LOGICAL ok_geost ! Initialisation vent geost. ou nul
72	LOGICAL ok_pv ! Polar Vortex
73	REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex
74
75	real zz,ran1
76	integer idum
77
78	REAL zdtvr
79
80	character(len=*),parameter :: modname="iniacademic"
81	character(len=80) :: abort_message
82
83	! Sanity check: verify that options selected by user are not incompatible
84	if ((iflag_phys==1).and. .not. read_start) then
85	write(lunout,*) trim(modname)," error: if read_start is set to ", &
86	" false then iflag_phys should not be 1"
87	write(lunout,*) "You most likely want an aquaplanet initialisation", &
88	" (iflag_phys >= 100)"
89	call abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.false.",1)
90	endif
91
92	!-----------------------------------------------------------------------
93	! 1. Initializations for Earth-like case
94	! --------------------------------------
95	!
96	! initialize planet radius, rotation rate,...
97	call conf_planete
98
99	time_0=0.
100	day_ref=1
101	annee_ref=0
102
103	im = iim
104	jm = jjm
105	day_ini = 1
106	dtvr = daysec/REAL(day_step)
107	zdtvr=dtvr
108	etot0 = 0.
109	ptot0 = 0.
110	ztot0 = 0.
111	stot0 = 0.
112	ang0 = 0.
113
114	if (llm == 1) then
115	! specific initializations for the shallow water case
116	kappa=1
117	endif
118
119	CALL iniconst
120	CALL inigeom
121	CALL inifilr
122
123	if (llm == 1) then
124	! initialize fields for the shallow water case, if required
125	if (.not.read_start) then
126	phis(ijb_u:ije_u)=0.
127	q(ijb_u:ije_u,1:llm,1:nqtot)=0
128	CALL sw_case_williamson91_6_loc(vcov,ucov,teta,masse,ps)
129	endif
130	endif
131
132	academic_case: if (iflag_phys >= 2) then
133	! initializations
134
135	! 1. local parameters
136	! by convention, winter is in the southern hemisphere
137	! Geostrophic wind or no wind?
138	ok_geost=.TRUE.
139	CALL getin('ok_geost',ok_geost)
140	! Constants for Newtonian relaxation and friction
141	k_f=1. !friction
142	CALL getin('k_j',k_f)
143	k_f=1./(daysec*k_f)
144	k_c_s=4. !cooling surface
145	CALL getin('k_c_s',k_c_s)
146	k_c_s=1./(daysec*k_c_s)
147	k_c_a=40. !cooling free atm
148	CALL getin('k_c_a',k_c_a)
149	k_c_a=1./(daysec*k_c_a)
150	! Constants for Teta equilibrium profile
151	teta0=315. ! mean Teta (S.H. 315K)
152	CALL getin('teta0',teta0)
153	ttp=200. ! Tropopause temperature (S.H. 200K)
154	CALL getin('ttp',ttp)
155	eps=0. ! Deviation to N-S symmetry(~0-20K)
156	CALL getin('eps',eps)
157	delt_y=60. ! Merid Temp. Gradient (S.H. 60K)
158	CALL getin('delt_y',delt_y)
159	delt_z=10. ! Vertical Gradient (S.H. 10K)
160	CALL getin('delt_z',delt_z)
161	! Polar vortex
162	ok_pv=.false.
163	CALL getin('ok_pv',ok_pv)
164	phi_pv=-50. ! Latitude of edge of vortex
165	CALL getin('phi_pv',phi_pv)
166	phi_pv=phi_pv*pi/180.
167	dphi_pv=5. ! Width of the edge
168	CALL getin('dphi_pv',dphi_pv)
169	dphi_pv=dphi_pv*pi/180.
170	gam_pv=4. ! -dT/dz vortex (in K/km)
171	CALL getin('gam_pv',gam_pv)
172
173	! 2. Initialize fields towards which to relax
174	! Friction
175	knewt_g=k_c_a
176	DO l=1,llm
177	zsig=presnivs(l)/preff
178	knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
179	kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
180	ENDDO
181	DO j=1,jjp1
182	clat4((j-1)iip1+1:jiip1)=cos(rlatu(j))**4
183	ENDDO
184
185	! Potential temperature
186	DO l=1,llm
187	zsig=presnivs(l)/preff
188	tetastrat=ttpzsig*(-kappa)
189	tetapv=tetastrat
190	IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
191	tetapv=tetastrat(zsig10.)*(kappacpp*gam_pv/1000./g)
192	ENDIF
193	DO j=1,jjp1
194	! Troposphere
195	ddsin=sin(rlatu(j))
196	tetajl(j,l)=teta0-delt_yddsinddsin+eps*ddsin &
197	-delt_z(1.-ddsinddsin)*log(zsig)
198	if (planet_type=="giant") then
199	tetajl(j,l)=teta0+(delt_y* &
200	((sin(rlatu(j)3.14159eps+0.0001))**2) &
201	/ ((rlatu(j)3.14159eps+0.0001)**2)) &
202	-delt_z*log(zsig)
203	endif
204	! Profil stratospherique isotherme (+vortex)
205	w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
206	tetastrat=tetastrat(1.-w_pv)+tetapvw_pv
207	tetajl(j,l)=MAX(tetajl(j,l),tetastrat)
208	ENDDO
209	ENDDO
210
211	! CALL writefield('theta_eq',tetajl)
212
213	do l=1,llm
214	do j=1,jjp1
215	do i=1,iip1
216	ij=(j-1)*iip1+i
217	tetarappel(ij,l)=tetajl(j,l)
218	enddo
219	enddo
220	enddo
221
222	! 3. Initialize fields (if necessary)
223	IF (.NOT. read_start) THEN
224	! allocate global fields:
225	! allocate(vcov_glo(ip1jm,llm))
226	allocate(ucov_glo(ip1jmp1,llm))
227	allocate(teta_glo(ip1jmp1,llm))
228	allocate(ps_glo(ip1jmp1))
229	allocate(masse_glo(ip1jmp1,llm))
230	allocate(phis_glo(ip1jmp1))
231
232	! surface pressure
233	if (iflag_phys>2) then
234	! specific value for CMIP5 aqua/terra planets
235	! "Specify the initial dry mass to be equivalent to
236	! a global mean surface pressure (101325 minus 245) Pa."
237	ps_glo(:)=101080.
238	else
239	! use reference surface pressure
240	ps_glo(:)=preff
241	endif
242
243	! ground geopotential
244	phis_glo(:)=0.
245
246	CALL pression ( ip1jmp1, ap, bp, ps_glo, p )
247	if (pressure_exner) then
248	CALL exner_hyb( ip1jmp1, ps_glo, p, pks, pk )
249	else
250	call exner_milieu(ip1jmp1,ps_glo,p,pks,pk)
251	endif
252	CALL massdair(p,masse_glo)
253
254	! bulk initialization of temperature
255	teta_glo(:,:)=tetarappel(:,:)
256
257	! geopotential
258	CALL geopot(ip1jmp1,teta_glo,pk,pks,phis_glo,phi)
259
260	! winds
261	if (ok_geost) then
262	call ugeostr(phi,ucov_glo)
263	else
264	ucov_glo(:,:)=0.
265	endif
266	vcov(ijb_v:ije_v,1:llm)=0.
267
268	! bulk initialization of tracers
269	if (planet_type=="earth") then
270	! Earth: first two tracers will be water
271	do i=1,nqtot
272	if (i == 1) q(ijb_u:ije_u,:,i)=1.e-10
273	if (i == 2) q(ijb_u:ije_u,:,i)=1.e-15
274	if (i.gt.2) q(ijb_u:ije_u,:,i)=0.
275	enddo
276	else
277	q(ijb_u:ije_u,:,:)=0
278	endif ! of if (planet_type=="earth")
279
280	! add random perturbation to temperature
281	idum = -1
282	zz = ran1(idum)
283	idum = 0
284	do l=1,llm
285	do ij=iip2,ip1jm
286	teta_glo(ij,l)=teta_glo(ij,l)(1.+0.005ran1(idum))
287	enddo
288	enddo
289
290	! maintain periodicity in longitude
291	do l=1,llm
292	do ij=1,ip1jmp1,iip1
293	teta_glo(ij+iim,l)=teta_glo(ij,l)
294	enddo
295	enddo
296
297	! copy data from global array to local array:
298	teta(ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:)
299	ucov(ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:)
300	! vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:)
301	masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:)
302	ps(ijb_u:ije_u)=ps_glo(ijb_u:ije_u)
303	phis(ijb_u:ije_u)=phis_glo(ijb_u:ije_u)
304
305	deallocate(teta_glo)
306	deallocate(ucov_glo)
307	! deallocate(vcov_glo)
308	deallocate(masse_glo)
309	deallocate(ps_glo)
310	deallocate(phis_glo)
311	ENDIF ! of IF (.NOT. read_start)
312	endif academic_case
313
314	END SUBROUTINE iniacademic_loc

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

Download in other formats:

Original Format