Context Navigation

iniacademic_loc.F90 @ 5134

Last change on this file since 5134 was 5084, checked in by Laurent Fairhead, 12 months ago
Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks AB, LF
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: 12.9 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, getKey, isoName
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 readTracFiles_mod, ONLY: addPhase
25	use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
26	use netcdf, only : NF90_CLOSE, NF90_GET_VAR
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,PARAMETER :: tnat1=.true.
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	! Earth: first two tracers will be water
314	do iq=1,nqtot
315	q(ijb_u:ije_u,:,iq)=0.
316	IF(tracers(iq)%name == addPhase('H2O', 'g')) q(ijb_u:ije_u,:,iq)=1.e-10
317	IF(tracers(iq)%name == addPhase('H2O', 'l')) q(ijb_u:ije_u,:,iq)=1.e-15
318
319	! CRisi: init des isotopes
320	! distill de Rayleigh très simplifiée
321	iName = tracers(iq)%iso_iName
322	if (niso <= 0 .OR. iName <= 0) CYCLE
323	iPhase = tracers(iq)%iso_iPhase
324	iqParent = tracers(iq)%iqParent
325	IF(tracers(iq)%iso_iZone == 0) THEN
326	if (tnat1) then
327	tnat=1.0
328	alpha_ideal=1.0
329	write(,) 'Attention dans iniacademic: alpha_ideal=1'
330	else
331	IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
332	CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
333	endif
334	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.)
335	ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
336	IF(tracers(iq)%iso_iZone == 1) THEN ! a verifier.
337	! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
338	! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
339	q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
340	else !IF(tracers(iq)%iso_iZone == 1) THEN
341	q(ijb_u:ije_u,:,iq) = 0.0
342	endif !IF(tracers(iq)%iso_iZone == 1) THEN
343	END IF !IF(tracers(iq)%iso_iZone == 0) THEN
344	enddo
345	else
346	q(ijb_u:ije_u,:,:)=0
347	endif ! of if (planet_type=="earth")
348
349	call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc')
350
351	! add random perturbation to temperature
352	idum = -1
353	zz = ran1(idum)
354	idum = 0
355	do l=1,llm
356	do ij=iip2,ip1jm
357	teta_glo(ij,l)=teta_glo(ij,l)(1.+tetanoiseran1(idum))
358	enddo
359	enddo
360
361	! maintain periodicity in longitude
362	do l=1,llm
363	do ij=1,ip1jmp1,iip1
364	teta_glo(ij+iim,l)=teta_glo(ij,l)
365	enddo
366	enddo
367
368	! copy data from global array to local array:
369	teta(ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:)
370	ucov(ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:)
371	! vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:)
372	masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:)
373	ps(ijb_u:ije_u)=ps_glo(ijb_u:ije_u)
374	phis(ijb_u:ije_u)=phis_glo(ijb_u:ije_u)
375
376	deallocate(teta_glo)
377	deallocate(ucov_glo)
378	! deallocate(vcov_glo)
379	deallocate(masse_glo)
380	deallocate(ps_glo)
381	deallocate(phis_glo)
382	ENDIF ! of IF (.NOT. read_start)
383	endif academic_case
384
385	END SUBROUTINE iniacademic_loc

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

Download in other formats:

Original Format