Context Navigation

iniacademic_loc.F90 @ 4249

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

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

Download in other formats:

Original Format