1 |
|
|
! |
2 |
|
|
! $Id: iniacademic.F90 4419 2023-02-05 20:32:11Z fhourdin $ |
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 |
|
|
|
83 |
|
|
character(len=*),parameter :: modname="iniacademic" |
84 |
|
|
character(len=80) :: abort_message |
85 |
|
|
|
86 |
|
|
! Sanity check: verify that options selected by user are not incompatible |
87 |
|
|
if ((iflag_phys==1).and. .not. read_start) then |
88 |
|
|
write(lunout,*) trim(modname)," error: if read_start is set to ", & |
89 |
|
|
" false then iflag_phys should not be 1" |
90 |
|
|
write(lunout,*) "You most likely want an aquaplanet initialisation", & |
91 |
|
|
" (iflag_phys >= 100)" |
92 |
|
|
call abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.false.",1) |
93 |
|
|
endif |
94 |
|
|
|
95 |
|
|
!----------------------------------------------------------------------- |
96 |
|
|
! 1. Initializations for Earth-like case |
97 |
|
|
! -------------------------------------- |
98 |
|
|
! |
99 |
|
|
! initialize planet radius, rotation rate,... |
100 |
|
|
call conf_planete |
101 |
|
|
|
102 |
|
|
time_0=0. |
103 |
|
|
day_ref=1 |
104 |
|
|
annee_ref=0 |
105 |
|
|
|
106 |
|
|
im = iim |
107 |
|
|
jm = jjm |
108 |
|
|
day_ini = 1 |
109 |
|
|
dtvr = daysec/REAL(day_step) |
110 |
|
|
zdtvr=dtvr |
111 |
|
|
etot0 = 0. |
112 |
|
|
ptot0 = 0. |
113 |
|
|
ztot0 = 0. |
114 |
|
|
stot0 = 0. |
115 |
|
|
ang0 = 0. |
116 |
|
|
|
117 |
|
|
if (llm == 1) then |
118 |
|
|
! specific initializations for the shallow water case |
119 |
|
|
kappa=1 |
120 |
|
|
endif |
121 |
|
|
|
122 |
|
|
CALL iniconst |
123 |
|
|
CALL inigeom |
124 |
|
|
CALL inifilr |
125 |
|
|
|
126 |
|
|
|
127 |
|
|
!------------------------------------------------------------------ |
128 |
|
|
! Initialize pressure and mass field if read_start=.false. |
129 |
|
|
!------------------------------------------------------------------ |
130 |
|
|
|
131 |
|
|
IF (.NOT. read_start) THEN |
132 |
|
|
|
133 |
|
|
!------------------------------------------------------------------ |
134 |
|
|
! Lecture eventuelle d'un fichier de relief interpollee sur la grille |
135 |
|
|
! du modele. |
136 |
|
|
! On suppose que le fichier relief_in.nc est stoké sur une grille |
137 |
|
|
! iim*jjp1 |
138 |
|
|
! Facile a créer à partir de la commande |
139 |
|
|
! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc |
140 |
|
|
!------------------------------------------------------------------ |
141 |
|
|
|
142 |
|
|
relief=0. |
143 |
|
|
ierr = NF90_OPEN ('relief_in.nc', NF90_NOWRITE,nid_relief) |
144 |
|
|
if (ierr.EQ.NF90_NOERR) THEN |
145 |
|
|
ierr=NF90_INQ_VARID(nid_relief,'RELIEF',varid) |
146 |
|
|
if (ierr==NF90_NOERR) THEN |
147 |
|
|
ierr=NF90_GET_VAR(nid_relief,varid,relief(1:iim,1:jjp1)) |
148 |
|
|
relief(iip1,:)=relief(1,:) |
149 |
|
|
else |
150 |
|
|
CALL abort_gcm ('iniacademic','variable RELIEF pas la',1) |
151 |
|
|
endif |
152 |
|
|
endif |
153 |
|
|
ierr = NF90_CLOSE (nid_relief) |
154 |
|
|
|
155 |
|
|
!------------------------------------------------------------------ |
156 |
|
|
! Initialisation du geopotentiel au sol et de la pression |
157 |
|
|
!------------------------------------------------------------------ |
158 |
|
|
|
159 |
|
|
print*,'relief=',minval(relief),maxval(relief),'g=',g |
160 |
|
|
do j=1,jjp1 |
161 |
|
|
do i=1,iip1 |
162 |
|
|
phis((j-1)*iip1+i)=g*relief(i,j) |
163 |
|
|
enddo |
164 |
|
|
enddo |
165 |
|
|
print*,'phis=',minval(phis),maxval(phis),'g=',g |
166 |
|
|
|
167 |
|
|
! ground geopotential |
168 |
|
|
!phis(:)=0. |
169 |
|
|
ps(:)=preff |
170 |
|
|
CALL pression ( ip1jmp1, ap, bp, ps, p ) |
171 |
|
|
|
172 |
|
|
if (pressure_exner) then |
173 |
|
|
CALL exner_hyb( ip1jmp1, ps, p, pks, pk) |
174 |
|
|
else |
175 |
|
|
call exner_milieu(ip1jmp1,ps,p,pks,pk) |
176 |
|
|
endif |
177 |
|
|
CALL massdair(p,masse) |
178 |
|
|
ENDIF |
179 |
|
|
|
180 |
|
|
if (llm == 1) then |
181 |
|
|
! initialize fields for the shallow water case, if required |
182 |
|
|
if (.not.read_start) then |
183 |
|
|
phis(:)=0. |
184 |
|
|
q(:,:,:)=0 |
185 |
|
|
CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) |
186 |
|
|
endif |
187 |
|
|
endif |
188 |
|
|
|
189 |
|
|
academic_case: if (iflag_phys >= 2) then |
190 |
|
|
! initializations |
191 |
|
|
|
192 |
|
|
! 1. local parameters |
193 |
|
|
! by convention, winter is in the southern hemisphere |
194 |
|
|
! Geostrophic wind or no wind? |
195 |
|
|
ok_geost=.TRUE. |
196 |
|
|
CALL getin('ok_geost',ok_geost) |
197 |
|
|
! Constants for Newtonian relaxation and friction |
198 |
|
|
k_f=1. !friction |
199 |
|
|
CALL getin('k_j',k_f) |
200 |
|
|
k_f=1./(daysec*k_f) |
201 |
|
|
k_c_s=4. !cooling surface |
202 |
|
|
CALL getin('k_c_s',k_c_s) |
203 |
|
|
k_c_s=1./(daysec*k_c_s) |
204 |
|
|
k_c_a=40. !cooling free atm |
205 |
|
|
CALL getin('k_c_a',k_c_a) |
206 |
|
|
k_c_a=1./(daysec*k_c_a) |
207 |
|
|
! Constants for Teta equilibrium profile |
208 |
|
|
teta0=315. ! mean Teta (S.H. 315K) |
209 |
|
|
CALL getin('teta0',teta0) |
210 |
|
|
ttp=200. ! Tropopause temperature (S.H. 200K) |
211 |
|
|
CALL getin('ttp',ttp) |
212 |
|
|
eps=0. ! Deviation to N-S symmetry(~0-20K) |
213 |
|
|
CALL getin('eps',eps) |
214 |
|
|
delt_y=60. ! Merid Temp. Gradient (S.H. 60K) |
215 |
|
|
CALL getin('delt_y',delt_y) |
216 |
|
|
delt_z=10. ! Vertical Gradient (S.H. 10K) |
217 |
|
|
CALL getin('delt_z',delt_z) |
218 |
|
|
! Polar vortex |
219 |
|
|
ok_pv=.false. |
220 |
|
|
CALL getin('ok_pv',ok_pv) |
221 |
|
|
phi_pv=-50. ! Latitude of edge of vortex |
222 |
|
|
CALL getin('phi_pv',phi_pv) |
223 |
|
|
phi_pv=phi_pv*pi/180. |
224 |
|
|
dphi_pv=5. ! Width of the edge |
225 |
|
|
CALL getin('dphi_pv',dphi_pv) |
226 |
|
|
dphi_pv=dphi_pv*pi/180. |
227 |
|
|
gam_pv=4. ! -dT/dz vortex (in K/km) |
228 |
|
|
CALL getin('gam_pv',gam_pv) |
229 |
|
|
tetanoise=0.005 |
230 |
|
|
CALL getin('tetanoise',tetanoise) |
231 |
|
|
|
232 |
|
|
! 2. Initialize fields towards which to relax |
233 |
|
|
! Friction |
234 |
|
|
knewt_g=k_c_a |
235 |
|
|
DO l=1,llm |
236 |
|
|
zsig=presnivs(l)/preff |
237 |
|
|
knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) |
238 |
|
|
kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) |
239 |
|
|
ENDDO |
240 |
|
|
DO j=1,jjp1 |
241 |
|
|
clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 |
242 |
|
|
ENDDO |
243 |
|
|
|
244 |
|
|
! Potential temperature |
245 |
|
|
DO l=1,llm |
246 |
|
|
zsig=presnivs(l)/preff |
247 |
|
|
tetastrat=ttp*zsig**(-kappa) |
248 |
|
|
tetapv=tetastrat |
249 |
|
|
IF ((ok_pv).AND.(zsig.LT.0.1)) THEN |
250 |
|
|
tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) |
251 |
|
|
ENDIF |
252 |
|
|
DO j=1,jjp1 |
253 |
|
|
! Troposphere |
254 |
|
|
ddsin=sin(rlatu(j)) |
255 |
|
|
tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin & |
256 |
|
|
-delt_z*(1.-ddsin*ddsin)*log(zsig) |
257 |
|
|
if (planet_type=="giant") then |
258 |
|
|
tetajl(j,l)=teta0+(delt_y* & |
259 |
|
|
((sin(rlatu(j)*3.14159*eps+0.0001))**2) & |
260 |
|
|
/ ((rlatu(j)*3.14159*eps+0.0001)**2)) & |
261 |
|
|
-delt_z*log(zsig) |
262 |
|
|
endif |
263 |
|
|
! Profil stratospherique isotherme (+vortex) |
264 |
|
|
w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. |
265 |
|
|
tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv |
266 |
|
|
tetajl(j,l)=MAX(tetajl(j,l),tetastrat) |
267 |
|
|
ENDDO |
268 |
|
|
ENDDO |
269 |
|
|
|
270 |
|
|
! CALL writefield('theta_eq',tetajl) |
271 |
|
|
|
272 |
|
|
do l=1,llm |
273 |
|
|
do j=1,jjp1 |
274 |
|
|
do i=1,iip1 |
275 |
|
|
ij=(j-1)*iip1+i |
276 |
|
|
tetarappel(ij,l)=tetajl(j,l) |
277 |
|
|
enddo |
278 |
|
|
enddo |
279 |
|
|
enddo |
280 |
|
|
|
281 |
|
|
! 3. Initialize fields (if necessary) |
282 |
|
|
IF (.NOT. read_start) THEN |
283 |
|
|
! bulk initialization of temperature |
284 |
|
|
IF (iflag_phys>10000) THEN |
285 |
|
|
! Particular case to impose a constant temperature T0=0.01*iflag_physx |
286 |
|
|
teta(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp) |
287 |
|
|
ELSE |
288 |
|
|
teta(:,:)=tetarappel(:,:) |
289 |
|
|
ENDIF |
290 |
|
|
! geopotential |
291 |
|
|
CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) |
292 |
|
|
|
293 |
|
|
DO l=1,llm |
294 |
|
|
print*,'presnivs,play,l',presnivs(l),(pk(1,l)/cpp)**(1./kappa)*preff |
295 |
|
|
!pks(ij) = (cpp/preff) * ps(ij) |
296 |
|
|
!pk(ij,1) = .5*pks(ij) |
297 |
|
|
! pk = cpp * (p/preff)^kappa |
298 |
|
|
ENDDO |
299 |
|
|
|
300 |
|
|
! winds |
301 |
|
|
if (ok_geost) then |
302 |
|
|
call ugeostr(phi,ucov) |
303 |
|
|
else |
304 |
|
|
ucov(:,:)=0. |
305 |
|
|
endif |
306 |
|
|
vcov(:,:)=0. |
307 |
|
|
|
308 |
|
|
! bulk initialization of tracers |
309 |
|
|
if (planet_type=="earth") then |
310 |
|
|
! Earth: first two tracers will be water |
311 |
|
|
do iq=1,nqtot |
312 |
|
|
q(:,:,iq)=0. |
313 |
|
|
IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10 |
314 |
|
|
IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15 |
315 |
|
|
|
316 |
|
|
! CRisi: init des isotopes |
317 |
|
|
! distill de Rayleigh très simplifiée |
318 |
|
|
iName = tracers(iq)%iso_iName |
319 |
|
|
if (niso <= 0 .OR. iName <= 0) CYCLE |
320 |
|
|
iPhase = tracers(iq)%iso_iPhase |
321 |
|
|
iqParent = tracers(iq)%iqParent |
322 |
|
|
IF(tracers(iq)%iso_iZone == 0) THEN |
323 |
|
|
IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & |
324 |
|
|
CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) |
325 |
|
|
q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.) |
326 |
|
|
ELSE |
327 |
|
|
q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase)) |
328 |
|
|
END IF |
329 |
|
|
enddo |
330 |
|
|
else |
331 |
|
|
q(:,:,:)=0 |
332 |
|
|
endif ! of if (planet_type=="earth") |
333 |
|
|
|
334 |
|
|
call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc') |
335 |
|
|
|
336 |
|
|
! add random perturbation to temperature |
337 |
|
|
idum = -1 |
338 |
|
|
zz = ran1(idum) |
339 |
|
|
idum = 0 |
340 |
|
|
do l=1,llm |
341 |
|
|
do ij=iip2,ip1jm |
342 |
|
|
teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum)) |
343 |
|
|
enddo |
344 |
|
|
enddo |
345 |
|
|
|
346 |
|
|
! maintain periodicity in longitude |
347 |
|
|
do l=1,llm |
348 |
|
|
do ij=1,ip1jmp1,iip1 |
349 |
|
|
teta(ij+iim,l)=teta(ij,l) |
350 |
|
|
enddo |
351 |
|
|
enddo |
352 |
|
|
|
353 |
|
|
ENDIF ! of IF (.NOT. read_start) |
354 |
|
|
endif academic_case |
355 |
|
|
|
356 |
|
|
END SUBROUTINE iniacademic |