Context Navigation

iniacademic_loc.F90 @ 5076

Last change on this file since 5076 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

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 lmdz_netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID,NF90_CLOSE, NF90_GET_VAR
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, iq, iName, iPhase, iqParent
70
71	integer :: nid_relief,varid,ierr
72	real, dimension(iip1,jjp1) :: relief
73
74
75	REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
76	REAL k_f,k_c_a,k_c_s ! Constantes de rappel
77	LOGICAL ok_geost ! Initialisation vent geost. ou nul
78	LOGICAL ok_pv ! Polar Vortex
79	REAL phi_pv,dphi_pv,gam_pv,tetanoise ! Constantes pour polar vortex
80
81	real zz,ran1
82	integer idum
83
84	REAL zdtvr, tnat, alpha_ideal
85	LOGICAL,PARAMETER :: tnat1=.true.
86
87	character(len=*),parameter :: modname="iniacademic"
88	character(len=80) :: abort_message
89
90	! Sanity check: verify that options selected by user are not incompatible
91	if ((iflag_phys==1).and. .not. read_start) then
92	write(lunout,*) trim(modname)," error: if read_start is set to ", &
93	" false then iflag_phys should not be 1"
94	write(lunout,*) "You most likely want an aquaplanet initialisation", &
95	" (iflag_phys >= 100)"
96	call abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.false.",1)
97	endif
98
99	!-----------------------------------------------------------------------
100	! 1. Initializations for Earth-like case
101	! --------------------------------------
102	!
103	! initialize planet radius, rotation rate,...
104	call conf_planete
105
106	time_0=0.
107	day_ref=1
108	! annee_ref=0
109
110	im = iim
111	jm = jjm
112	day_ini = 1
113	dtvr = daysec/REAL(day_step)
114	zdtvr=dtvr
115	etot0 = 0.
116	ptot0 = 0.
117	ztot0 = 0.
118	stot0 = 0.
119	ang0 = 0.
120
121	if (llm == 1) then
122	! specific initializations for the shallow water case
123	kappa=1
124	endif
125
126	CALL iniconst
127	CALL inigeom
128	CALL inifilr
129
130	! Initialize pressure and mass field if read_start=.false.
131	IF (.NOT. read_start) THEN
132	! allocate global fields:
133	! allocate(vcov_glo(ip1jm,llm))
134
135	allocate(ucov_glo(ip1jmp1,llm))
136	allocate(teta_glo(ip1jmp1,llm))
137	allocate(ps_glo(ip1jmp1))
138	allocate(masse_glo(ip1jmp1,llm))
139	allocate(phis_glo(ip1jmp1))
140
141	! surface pressure
142	ps_glo(:)=preff
143
144	!------------------------------------------------------------------
145	! Lecture eventuelle d'un fichier de relief interpollee sur la grille
146	! du modele.
147	! On suppose que le fichier relief_in.nc est stoké sur une grille
148	! iim*jjp1
149	! Facile a créer à partir de la commande
150	! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
151	!------------------------------------------------------------------
152
153	relief=0.
154	ierr = NF90_OPEN ('relief_in.nc', NF90_NOWRITE,nid_relief)
155	if (ierr==NF90_NOERR) THEN
156	ierr=NF90_INQ_VARID(nid_relief,'RELIEF',varid)
157	if (ierr==NF90_NOERR) THEN
158	ierr=NF90_GET_VAR(nid_relief,varid,relief(1:iim,1:jjp1))
159	relief(iip1,:)=relief(1,:)
160	else
161	CALL abort_gcm ('iniacademic','variable RELIEF pas la',1)
162	endif
163	endif
164	ierr = NF90_CLOSE (nid_relief)
165
166
167	!------------------------------------------------------------------
168	! Initialisation du geopotentiel au sol et de la pression
169	!------------------------------------------------------------------
170
171	print*,'relief=',minval(relief),maxval(relief),'g=',g
172	do j=1,jjp1
173	do i=1,iip1
174	phis_glo((j-1)iip1+i)=grelief(i,j)
175	enddo
176	enddo
177	print*,'phis=',minval(phis),maxval(phis),'g=',g
178
179	CALL pression ( ip1jmp1, ap, bp, ps_glo, p )
180	if (pressure_exner) then
181	CALL exner_hyb( ip1jmp1, ps_glo, p, pks, pk )
182	else
183	call exner_milieu(ip1jmp1,ps_glo,p,pks,pk)
184	endif
185	CALL massdair(p,masse_glo)
186	ENDIF
187
188	if (llm == 1) then
189	! initialize fields for the shallow water case, if required
190	if (.not.read_start) then
191	phis(ijb_u:ije_u)=0.
192	q(ijb_u:ije_u,1:llm,1:nqtot)=0
193	CALL sw_case_williamson91_6_loc(vcov,ucov,teta,masse,ps)
194	endif
195	endif
196
197	academic_case: if (iflag_phys >= 2) then
198	! initializations
199
200	! 1. local parameters
201	! by convention, winter is in the southern hemisphere
202	! Geostrophic wind or no wind?
203	ok_geost=.TRUE.
204	CALL getin('ok_geost',ok_geost)
205	! Constants for Newtonian relaxation and friction
206	k_f=1. !friction
207	CALL getin('k_j',k_f)
208	k_f=1./(daysec*k_f)
209	k_c_s=4. !cooling surface
210	CALL getin('k_c_s',k_c_s)
211	k_c_s=1./(daysec*k_c_s)
212	k_c_a=40. !cooling free atm
213	CALL getin('k_c_a',k_c_a)
214	k_c_a=1./(daysec*k_c_a)
215	! Constants for Teta equilibrium profile
216	teta0=315. ! mean Teta (S.H. 315K)
217	CALL getin('teta0',teta0)
218	ttp=200. ! Tropopause temperature (S.H. 200K)
219	CALL getin('ttp',ttp)
220	eps=0. ! Deviation to N-S symmetry(~0-20K)
221	CALL getin('eps',eps)
222	delt_y=60. ! Merid Temp. Gradient (S.H. 60K)
223	CALL getin('delt_y',delt_y)
224	delt_z=10. ! Vertical Gradient (S.H. 10K)
225	CALL getin('delt_z',delt_z)
226	! Polar vortex
227	ok_pv=.false.
228	CALL getin('ok_pv',ok_pv)
229	phi_pv=-50. ! Latitude of edge of vortex
230	CALL getin('phi_pv',phi_pv)
231	phi_pv=phi_pv*pi/180.
232	dphi_pv=5. ! Width of the edge
233	CALL getin('dphi_pv',dphi_pv)
234	dphi_pv=dphi_pv*pi/180.
235	gam_pv=4. ! -dT/dz vortex (in K/km)
236	CALL getin('gam_pv',gam_pv)
237	tetanoise=0.005
238	CALL getin('tetanoise',tetanoise)
239
240	! 2. Initialize fields towards which to relax
241	! Friction
242	knewt_g=k_c_a
243	DO l=1,llm
244	zsig=presnivs(l)/preff
245	knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
246	kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
247	ENDDO
248	DO j=1,jjp1
249	clat4((j-1)iip1+1:jiip1)=cos(rlatu(j))**4
250	ENDDO
251
252	! Potential temperature
253	DO l=1,llm
254	zsig=presnivs(l)/preff
255	tetastrat=ttpzsig*(-kappa)
256	tetapv=tetastrat
257	IF ((ok_pv).AND.(zsig<0.1)) THEN
258	tetapv=tetastrat(zsig10.)*(kappacpp*gam_pv/1000./g)
259	ENDIF
260	DO j=1,jjp1
261	! Troposphere
262	ddsin=sin(rlatu(j))
263	tetajl(j,l)=teta0-delt_yddsinddsin+eps*ddsin &
264	-delt_z(1.-ddsinddsin)*log(zsig)
265	if (planet_type=="giant") then
266	tetajl(j,l)=teta0+(delt_y* &
267	((sin(rlatu(j)3.14159eps+0.0001))**2) &
268	/ ((rlatu(j)3.14159eps+0.0001)**2)) &
269	-delt_z*log(zsig)
270	endif
271	! Profil stratospherique isotherme (+vortex)
272	w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
273	tetastrat=tetastrat(1.-w_pv)+tetapvw_pv
274	tetajl(j,l)=MAX(tetajl(j,l),tetastrat)
275	ENDDO
276	ENDDO
277
278	! CALL writefield('theta_eq',tetajl)
279
280	do l=1,llm
281	do j=1,jjp1
282	do i=1,iip1
283	ij=(j-1)*iip1+i
284	tetarappel(ij,l)=tetajl(j,l)
285	enddo
286	enddo
287	enddo
288
289	! 3. Initialize fields (if necessary)
290	IF (.NOT. read_start) THEN
291	! bulk initialization of temperature
292	IF (iflag_phys>10000) THEN
293	! Particular case to impose a constant temperature T0=0.01*iflag_phys
294	teta_glo(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp)
295	ELSE
296	teta_glo(:,:)=tetarappel(:,:)
297	ENDIF
298	! geopotential
299	CALL geopot(ip1jmp1,teta_glo,pk,pks,phis_glo,phi)
300
301	! winds
302	if (ok_geost) then
303	call ugeostr(phi,ucov_glo)
304	else
305	ucov_glo(:,:)=0.
306	endif
307	vcov(ijb_v:ije_v,1:llm)=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(ijb_u:ije_u,:,iq)=0.
314	IF(tracers(iq)%name == addPhase('H2O', 'g')) q(ijb_u:ije_u,:,iq)=1.e-10
315	IF(tracers(iq)%name == addPhase('H2O', 'l')) q(ijb_u:ije_u,:,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(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)tnat(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.)
333	ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
334	IF(tracers(iq)%iso_iZone == 1) THEN ! a verifier.
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(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
338	else !IF(tracers(iq)%iso_iZone == 1) THEN
339	q(ijb_u:ije_u,:,iq) = 0.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(ijb_u:ije_u,:,:)=0
345	endif ! of if (planet_type=="earth")
346
347	call check_isotopes(q,ijb_u,ije_u,'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_glo(ij,l)=teta_glo(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_glo(ij+iim,l)=teta_glo(ij,l)
363	enddo
364	enddo
365
366	! copy data from global array to local array:
367	teta(ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:)
368	ucov(ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:)
369	! vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:)
370	masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:)
371	ps(ijb_u:ije_u)=ps_glo(ijb_u:ije_u)
372	phis(ijb_u:ije_u)=phis_glo(ijb_u:ije_u)
373
374	deallocate(teta_glo)
375	deallocate(ucov_glo)
376	! deallocate(vcov_glo)
377	deallocate(masse_glo)
378	deallocate(ps_glo)
379	deallocate(phis_glo)
380	ENDIF ! of IF (.NOT. read_start)
381	endif academic_case
382
383	END SUBROUTINE iniacademic_loc

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

Download in other formats:

Original Format