Context Navigation

iniacademic.F90 @ 5073

Last change on this file since 5073 was 4984, checked in by crisi, 6 months ago
plenty of files that I forgot to commit last time.
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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 12.0 KB

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

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

Download in other formats:

Original Format