Context Navigation

iniacademic_loc.F90 @ 5227

Last change on this file since 5227 was 5214, checked in by dcugnet, 3 months ago
The fortran parameters file "iso_params_mod.F90" is introduced so that "tnat" and "alpha_ideal" are defined in a single place but used in several. The "getKey" routine is only used in "infotrac" and "infotrac_phy" routines, but could be used outside this scope to get tracers parameters (read from "tracer.def") or isotopic parameters (read from "isotopes_params.def" - disabled for now).
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: 13.4 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 infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase
8	USE control_mod, ONLY: day_step,planet_type
9	use exner_hyb_m, only: exner_hyb
10	use exner_milieu_m, only: exner_milieu
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	USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm
20	USE logic_mod, ONLY: iflag_phys, read_start
21	USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
22	USE temps_mod, ONLY: annee_ref, day_ini, day_ref
23	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
24	use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
25	use netcdf, only : NF90_CLOSE, NF90_GET_VAR
26	USE iso_params_mod ! tnat_* and alpha_ideal_*
27
28
29	! Author: Frederic Hourdin original: 15/01/93
30	! The forcing defined here is from Held and Suarez, 1994, Bulletin
31	! of the American Meteorological Society, 75, 1825.
32
33	IMPLICIT NONE
34
35	! Declararations:
36	! ---------------
37
38	include "dimensions.h"
39	include "paramet.h"
40	include "comgeom.h"
41	include "academic.h"
42	include "iniprint.h"
43
44	! Arguments:
45	! ----------
46
47	REAL,INTENT(OUT) :: time_0
48
49	! fields
50	REAL,INTENT(OUT) :: vcov(ijb_v:ije_v,llm) ! meridional covariant wind
51	REAL,INTENT(OUT) :: ucov(ijb_u:ije_u,llm) ! zonal covariant wind
52	REAL,INTENT(OUT) :: teta(ijb_u:ije_u,llm) ! potential temperature (K)
53	REAL,INTENT(OUT) :: q(ijb_u:ije_u,llm,nqtot) ! advected tracers (.../kg_of_air)
54	REAL,INTENT(OUT) :: ps(ijb_u:ije_u) ! surface pressure (Pa)
55	REAL,INTENT(OUT) :: masse(ijb_u:ije_u,llm) ! air mass in grid cell (kg)
56	REAL,INTENT(OUT) :: phis(ijb_u:ije_u) ! surface geopotential
57
58	! Local:
59	! ------
60
61	REAL,ALLOCATABLE :: vcov_glo(:,:),ucov_glo(:,:),teta_glo(:,:)
62	REAL,ALLOCATABLE :: q_glo(:,:),masse_glo(:,:),ps_glo(:)
63	REAL,ALLOCATABLE :: phis_glo(:)
64	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
65	REAL pks(ip1jmp1) ! exner au sol
66	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
67	REAL phi(ip1jmp1,llm) ! geopotentiel
68	REAL ddsin,zsig,tetapv,w_pv ! variables auxiliaires
69	real tetastrat ! potential temperature in the stratosphere, in K
70	real tetajl(jjp1,llm)
71	INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
72
73	integer :: nid_relief,varid,ierr
74	real, dimension(iip1,jjp1) :: relief
75
76
77	REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
78	REAL k_f,k_c_a,k_c_s ! Constantes de rappel
79	LOGICAL ok_geost ! Initialisation vent geost. ou nul
80	LOGICAL ok_pv ! Polar Vortex
81	REAL phi_pv,dphi_pv,gam_pv,tetanoise ! Constantes pour polar vortex
82
83	real zz,ran1
84	integer idum
85
86	REAL zdtvr, tnat, alpha_ideal
87	LOGICAL :: ltnat1
88
89	character(len=*),parameter :: modname="iniacademic"
90	character(len=80) :: abort_message
91
92	! Sanity check: verify that options selected by user are not incompatible
93	if ((iflag_phys==1).and. .not. read_start) then
94	write(lunout,*) trim(modname)," error: if read_start is set to ", &
95	" false then iflag_phys should not be 1"
96	write(lunout,*) "You most likely want an aquaplanet initialisation", &
97	" (iflag_phys >= 100)"
98	call abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.false.",1)
99	endif
100
101	!-----------------------------------------------------------------------
102	! 1. Initializations for Earth-like case
103	! --------------------------------------
104	!
105	! initialize planet radius, rotation rate,...
106	call conf_planete
107
108	time_0=0.
109	day_ref=1
110	! annee_ref=0
111
112	im = iim
113	jm = jjm
114	day_ini = 1
115	dtvr = daysec/REAL(day_step)
116	zdtvr=dtvr
117	etot0 = 0.
118	ptot0 = 0.
119	ztot0 = 0.
120	stot0 = 0.
121	ang0 = 0.
122
123	if (llm == 1) then
124	! specific initializations for the shallow water case
125	kappa=1
126	endif
127
128	CALL iniconst
129	CALL inigeom
130	CALL inifilr
131
132	! Initialize pressure and mass field if read_start=.false.
133	IF (.NOT. read_start) THEN
134	! allocate global fields:
135	! allocate(vcov_glo(ip1jm,llm))
136
137	allocate(ucov_glo(ip1jmp1,llm))
138	allocate(teta_glo(ip1jmp1,llm))
139	allocate(ps_glo(ip1jmp1))
140	allocate(masse_glo(ip1jmp1,llm))
141	allocate(phis_glo(ip1jmp1))
142
143	! surface pressure
144	ps_glo(:)=preff
145
146	!------------------------------------------------------------------
147	! Lecture eventuelle d'un fichier de relief interpollee sur la grille
148	! du modele.
149	! On suppose que le fichier relief_in.nc est stoké sur une grille
150	! iim*jjp1
151	! Facile a créer à partir de la commande
152	! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
153	!------------------------------------------------------------------
154
155	relief=0.
156	ierr = NF90_OPEN ('relief_in.nc', NF90_NOWRITE,nid_relief)
157	if (ierr.EQ.NF90_NOERR) THEN
158	ierr=NF90_INQ_VARID(nid_relief,'RELIEF',varid)
159	if (ierr==NF90_NOERR) THEN
160	ierr=NF90_GET_VAR(nid_relief,varid,relief(1:iim,1:jjp1))
161	relief(iip1,:)=relief(1,:)
162	else
163	CALL abort_gcm ('iniacademic','variable RELIEF pas la',1)
164	endif
165	endif
166	ierr = NF90_CLOSE (nid_relief)
167
168
169	!------------------------------------------------------------------
170	! Initialisation du geopotentiel au sol et de la pression
171	!------------------------------------------------------------------
172
173	print*,'relief=',minval(relief),maxval(relief),'g=',g
174	do j=1,jjp1
175	do i=1,iip1
176	phis_glo((j-1)iip1+i)=grelief(i,j)
177	enddo
178	enddo
179	print*,'phis=',minval(phis),maxval(phis),'g=',g
180
181	CALL pression ( ip1jmp1, ap, bp, ps_glo, p )
182	if (pressure_exner) then
183	CALL exner_hyb( ip1jmp1, ps_glo, p, pks, pk )
184	else
185	call exner_milieu(ip1jmp1,ps_glo,p,pks,pk)
186	endif
187	CALL massdair(p,masse_glo)
188	ENDIF
189
190	if (llm == 1) then
191	! initialize fields for the shallow water case, if required
192	if (.not.read_start) then
193	phis(ijb_u:ije_u)=0.
194	q(ijb_u:ije_u,1:llm,1:nqtot)=0
195	CALL sw_case_williamson91_6_loc(vcov,ucov,teta,masse,ps)
196	endif
197	endif
198
199	academic_case: if (iflag_phys >= 2) then
200	! initializations
201
202	! 1. local parameters
203	! by convention, winter is in the southern hemisphere
204	! Geostrophic wind or no wind?
205	ok_geost=.TRUE.
206	CALL getin('ok_geost',ok_geost)
207	! Constants for Newtonian relaxation and friction
208	k_f=1. !friction
209	CALL getin('k_j',k_f)
210	k_f=1./(daysec*k_f)
211	k_c_s=4. !cooling surface
212	CALL getin('k_c_s',k_c_s)
213	k_c_s=1./(daysec*k_c_s)
214	k_c_a=40. !cooling free atm
215	CALL getin('k_c_a',k_c_a)
216	k_c_a=1./(daysec*k_c_a)
217	! Constants for Teta equilibrium profile
218	teta0=315. ! mean Teta (S.H. 315K)
219	CALL getin('teta0',teta0)
220	ttp=200. ! Tropopause temperature (S.H. 200K)
221	CALL getin('ttp',ttp)
222	eps=0. ! Deviation to N-S symmetry(~0-20K)
223	CALL getin('eps',eps)
224	delt_y=60. ! Merid Temp. Gradient (S.H. 60K)
225	CALL getin('delt_y',delt_y)
226	delt_z=10. ! Vertical Gradient (S.H. 10K)
227	CALL getin('delt_z',delt_z)
228	! Polar vortex
229	ok_pv=.false.
230	CALL getin('ok_pv',ok_pv)
231	phi_pv=-50. ! Latitude of edge of vortex
232	CALL getin('phi_pv',phi_pv)
233	phi_pv=phi_pv*pi/180.
234	dphi_pv=5. ! Width of the edge
235	CALL getin('dphi_pv',dphi_pv)
236	dphi_pv=dphi_pv*pi/180.
237	gam_pv=4. ! -dT/dz vortex (in K/km)
238	CALL getin('gam_pv',gam_pv)
239	tetanoise=0.005
240	CALL getin('tetanoise',tetanoise)
241
242	! 2. Initialize fields towards which to relax
243	! Friction
244	knewt_g=k_c_a
245	DO l=1,llm
246	zsig=presnivs(l)/preff
247	knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
248	kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
249	ENDDO
250	DO j=1,jjp1
251	clat4((j-1)iip1+1:jiip1)=cos(rlatu(j))**4
252	ENDDO
253
254	! Potential temperature
255	DO l=1,llm
256	zsig=presnivs(l)/preff
257	tetastrat=ttpzsig*(-kappa)
258	tetapv=tetastrat
259	IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
260	tetapv=tetastrat(zsig10.)*(kappacpp*gam_pv/1000./g)
261	ENDIF
262	DO j=1,jjp1
263	! Troposphere
264	ddsin=sin(rlatu(j))
265	tetajl(j,l)=teta0-delt_yddsinddsin+eps*ddsin &
266	-delt_z(1.-ddsinddsin)*log(zsig)
267	if (planet_type=="giant") then
268	tetajl(j,l)=teta0+(delt_y* &
269	((sin(rlatu(j)3.14159eps+0.0001))**2) &
270	/ ((rlatu(j)3.14159eps+0.0001)**2)) &
271	-delt_z*log(zsig)
272	endif
273	! Profil stratospherique isotherme (+vortex)
274	w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
275	tetastrat=tetastrat(1.-w_pv)+tetapvw_pv
276	tetajl(j,l)=MAX(tetajl(j,l),tetastrat)
277	ENDDO
278	ENDDO
279
280	! CALL writefield('theta_eq',tetajl)
281
282	do l=1,llm
283	do j=1,jjp1
284	do i=1,iip1
285	ij=(j-1)*iip1+i
286	tetarappel(ij,l)=tetajl(j,l)
287	enddo
288	enddo
289	enddo
290
291	! 3. Initialize fields (if necessary)
292	IF (.NOT. read_start) THEN
293	! bulk initialization of temperature
294	IF (iflag_phys>10000) THEN
295	! Particular case to impose a constant temperature T0=0.01*iflag_phys
296	teta_glo(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp)
297	ELSE
298	teta_glo(:,:)=tetarappel(:,:)
299	ENDIF
300	! geopotential
301	CALL geopot(ip1jmp1,teta_glo,pk,pks,phis_glo,phi)
302
303	! winds
304	if (ok_geost) then
305	call ugeostr(phi,ucov_glo)
306	else
307	ucov_glo(:,:)=0.
308	endif
309	vcov(ijb_v:ije_v,1:llm)=0.
310
311	! bulk initialization of tracers
312	if (planet_type=="earth") then
313	ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
314	! Earth: first two tracers will be water
315	do iq=1,nqtot
316	q(ijb_u:ije_u,:,iq)=0.
317	IF(tracers(iq)%name == addPhase('H2O', 'g')) q(ijb_u:ije_u,:,iq)=1.e-10
318	IF(tracers(iq)%name == addPhase('H2O', 'l')) q(ijb_u:ije_u,:,iq)=1.e-15
319
320	! CRisi: init des isotopes
321	! distill de Rayleigh très simplifiée
322	iName = tracers(iq)%iso_iName
323	if (niso <= 0 .OR. iName <= 0) CYCLE
324	iPhase = tracers(iq)%iso_iPhase
325	iqParent = tracers(iq)%iqParent
326	IF(tracers(iq)%iso_iZone == 0) THEN
327	IF(ltnat1) THEN
328	tnat = 1.0
329	alpha_ideal = 1.0
330	WRITE(lunout, *) 'In '//TRIM(modname)//': !!! Beware: alpha_ideal put to 1 !!!'
331	ELSE
332	SELECT CASE(isoName(iName))
333	CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O
334	CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O
335	CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O
336	CASE('HDO'); tnat = tnat_HDO; alpha_ideal = alpha_ideal_HDO
337	CASE('HTO'); tnat = tnat_HTO; alpha_ideal = alpha_ideal_HTO
338	CASE DEFAULT
339	CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1)
340	END SELECT
341	END IF
342	q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)tnat(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.)
343	ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
344	IF(tracers(iq)%iso_iZone == 1) THEN ! a verifier.
345	! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
346	! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
347	q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
348	else !IF(tracers(iq)%iso_iZone == 1) THEN
349	q(ijb_u:ije_u,:,iq) = 0.0
350	endif !IF(tracers(iq)%iso_iZone == 1) THEN
351	END IF !IF(tracers(iq)%iso_iZone == 0) THEN
352	enddo
353	else
354	q(ijb_u:ije_u,:,:)=0
355	endif ! of if (planet_type=="earth")
356
357	call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc')
358
359	! add random perturbation to temperature
360	idum = -1
361	zz = ran1(idum)
362	idum = 0
363	do l=1,llm
364	do ij=iip2,ip1jm
365	teta_glo(ij,l)=teta_glo(ij,l)(1.+tetanoiseran1(idum))
366	enddo
367	enddo
368
369	! maintain periodicity in longitude
370	do l=1,llm
371	do ij=1,ip1jmp1,iip1
372	teta_glo(ij+iim,l)=teta_glo(ij,l)
373	enddo
374	enddo
375
376	! copy data from global array to local array:
377	teta(ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:)
378	ucov(ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:)
379	! vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:)
380	masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:)
381	ps(ijb_u:ije_u)=ps_glo(ijb_u:ije_u)
382	phis(ijb_u:ije_u)=phis_glo(ijb_u:ije_u)
383
384	deallocate(teta_glo)
385	deallocate(ucov_glo)
386	! deallocate(vcov_glo)
387	deallocate(masse_glo)
388	deallocate(ps_glo)
389	deallocate(phis_glo)
390	ENDIF ! of IF (.NOT. read_start)
391	endif academic_case
392
393	END SUBROUTINE iniacademic_loc

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

Download in other formats:

Original Format