Context Navigation

iniacademic.F90 @ 5080

Last change on this file since 5080 was 5075, checked in by abarral, 4 months ago

[continued & end] replace netcdf by lmdz_netcdf.F90 wrapper
"use netcdf" is now only used in lmdz_netcdf.F90 (except ecrad and obsolete/)
<include "netcdf.inc"> is now likewise only used in lmdz_netcdf.F90.

systematically specify explicitely <USE lmdz_netcdf, ONLY:> (probably left some missing, to correct later on)

Further replacement of nf_put_* by nf90_put_* (same for _get_)

[minor] replace deprecated boolean operators along the way

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

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

Download in other formats:

Original Format