My Project
 All Classes Files Functions Variables Macros
iniacademic_loc.F90
Go to the documentation of this file.
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
7  USE infotrac, ONLY : nqtot
8  USE control_mod, ONLY: day_step,planet_type
9  USE parallel
10 #ifdef CPP_IOIPSL
11  USE ioipsl
12 #else
13  ! if not using IOIPSL, we still need to use (a local version of) getin
14  USE ioipsl_getincom
15 #endif
16  USE write_field
17 
18  ! Author: Frederic Hourdin original: 15/01/93
19  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
20  ! of the American Meteorological Society, 75, 1825.
21 
22  IMPLICIT NONE
23 
24  ! Declararations:
25  ! ---------------
26 
27  include "dimensions.h"
28  include "paramet.h"
29  include "comvert.h"
30  include "comconst.h"
31  include "comgeom.h"
32  include "academic.h"
33  include "ener.h"
34  include "temps.h"
35  include "iniprint.h"
36  include "logic.h"
37 
38  ! Arguments:
39  ! ----------
40 
41  real time_0
42 
43  ! variables dynamiques
44  REAL vcov(ijb_v:ije_v,llm),ucov(ijb_u:ije_u,llm) ! vents covariants
45  REAL teta(ijb_u:ije_u,llm) ! temperature potentielle
46  REAL q(ijb_u:ije_u,llm,nqtot) ! champs advectes
47  REAL ps(ijb_u:ije_u) ! pression au sol
48  REAL masse(ijb_u:ije_u,llm) ! masse d'air
49  REAL phis(ijb_u:ije_u) ! geopotentiel au sol
50 
51  ! Local:
52  ! ------
53 
54  REAL,ALLOCATABLE :: vcov_glo(:,:),ucov_glo(:,:),teta_glo(:,:)
55  REAL,ALLOCATABLE :: q_glo(:,:),masse_glo(:,:),ps_glo(:)
56  REAL,ALLOCATABLE :: phis_glo(:)
57  REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
58  REAL pks(ip1jmp1) ! exner au sol
59  REAL pk(ip1jmp1,llm) ! exner au milieu des couches
60  REAL pkf(ip1jmp1,llm) ! exner filt.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
66 
67  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
68  REAL k_f,k_c_a,k_c_s ! Constantes de rappel
69  LOGICAL ok_geost ! Initialisation vent geost. ou nul
70  LOGICAL ok_pv ! Polar Vortex
71  REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex
72 
73  real zz,ran1
74  integer idum
75 
76  REAL zdtvr
77  real,allocatable :: alpha(:,:),beta(:,:)
78 
79  character(len=*),parameter :: modname="iniacademic"
80  character(len=80) :: abort_message
81 
82  !-----------------------------------------------------------------------
83  ! 1. Initializations for Earth-like case
84  ! --------------------------------------
85  !
86  ! initialize planet radius, rotation rate,...
87  call conf_planete
88 
89  time_0=0.
90  day_ref=1
91  annee_ref=0
92 
93  im = iim
94  jm = jjm
95  day_ini = 1
96  dtvr = daysec/REAL(day_step)
97  zdtvr=dtvr
98  etot0 = 0.
99  ptot0 = 0.
100  ztot0 = 0.
101  stot0 = 0.
102  ang0 = 0.
103 
104  if (llm == 1) then
105  ! specific initializations for the shallow water case
106  kappa=1
107  endif
108 
109  CALL iniconst
110  CALL inigeom
111  CALL inifilr
112 
113  if (llm == 1) then
114  ! initialize fields for the shallow water case, if required
115  if (.not.read_start) then
116  phis(ijb_u:ije_u)=0.
117  q(ijb_u:ije_u,1:llm,1:nqtot)=0
118  CALL sw_case_williamson91_6_loc(vcov,ucov,teta,masse,ps)
119  endif
120  endif
121 
122  academic_case: if (iflag_phys >= 2) then
123  ! initializations
124 
125  ! 1. local parameters
126  ! by convention, winter is in the southern hemisphere
127  ! Geostrophic wind or no wind?
128  ok_geost=.true.
129  CALL getin('ok_geost',ok_geost)
130  ! Constants for Newtonian relaxation and friction
131  k_f=1. !friction
132  CALL getin('k_j',k_f)
133  k_f=1./(daysec*k_f)
134  k_c_s=4. !cooling surface
135  CALL getin('k_c_s',k_c_s)
136  k_c_s=1./(daysec*k_c_s)
137  k_c_a=40. !cooling free atm
138  CALL getin('k_c_a',k_c_a)
139  k_c_a=1./(daysec*k_c_a)
140  ! Constants for Teta equilibrium profile
141  teta0=315. ! mean Teta (S.H. 315K)
142  CALL getin('teta0',teta0)
143  ttp=200. ! Tropopause temperature (S.H. 200K)
144  CALL getin('ttp',ttp)
145  eps=0. ! Deviation to N-S symmetry(~0-20K)
146  CALL getin('eps',eps)
147  delt_y=60. ! Merid Temp. Gradient (S.H. 60K)
148  CALL getin('delt_y',delt_y)
149  delt_z=10. ! Vertical Gradient (S.H. 10K)
150  CALL getin('delt_z',delt_z)
151  ! Polar vortex
152  ok_pv=.false.
153  CALL getin('ok_pv',ok_pv)
154  phi_pv=-50. ! Latitude of edge of vortex
155  CALL getin('phi_pv',phi_pv)
156  phi_pv=phi_pv*pi/180.
157  dphi_pv=5. ! Width of the edge
158  CALL getin('dphi_pv',dphi_pv)
159  dphi_pv=dphi_pv*pi/180.
160  gam_pv=4. ! -dT/dz vortex (in K/km)
161  CALL getin('gam_pv',gam_pv)
162 
163  ! 2. Initialize fields towards which to relax
164  ! Friction
165  knewt_g=k_c_a
166  DO l=1,llm
167  zsig=presnivs(l)/preff
168  knewt_t(l)=(k_c_s-k_c_a)*max(0.,(zsig-0.7)/0.3)
169  kfrict(l)=k_f*max(0.,(zsig-0.7)/0.3)
170  ENDDO
171  DO j=1,jjp1
172  clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
173  ENDDO
174 
175  ! Potential temperature
176  DO l=1,llm
177  zsig=presnivs(l)/preff
178  tetastrat=ttp*zsig**(-kappa)
179  tetapv=tetastrat
180  IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
181  tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
182  ENDIF
183  DO j=1,jjp1
184  ! Troposphere
185  ddsin=sin(rlatu(j))
186  tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
187  -delt_z*(1.-ddsin*ddsin)*log(zsig)
188  if (planet_type=="giant") then
189  tetajl(j,l)=teta0+(delt_y* &
190  ((sin(rlatu(j)*3.14159*eps+0.0001))**2) &
191  / ((rlatu(j)*3.14159*eps+0.0001)**2)) &
192  -delt_z*log(zsig)
193  endif
194  ! Profil stratospherique isotherme (+vortex)
195  w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
196  tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv
197  tetajl(j,l)=max(tetajl(j,l),tetastrat)
198  ENDDO
199  ENDDO
200 
201  ! CALL writefield('theta_eq',tetajl)
202 
203  do l=1,llm
204  do j=1,jjp1
205  do i=1,iip1
206  ij=(j-1)*iip1+i
207  tetarappel(ij,l)=tetajl(j,l)
208  enddo
209  enddo
210  enddo
211 
212  ! 3. Initialize fields (if necessary)
213  IF (.NOT. read_start) THEN
214  ! allocate global fields:
215 ! allocate(vcov_glo(ip1jm,llm))
216  allocate(ucov_glo(ip1jmp1,llm))
217  allocate(teta_glo(ip1jmp1,llm))
218  allocate(ps_glo(ip1jmp1))
219  allocate(masse_glo(ip1jmp1,llm))
220  allocate(phis_glo(ip1jmp1))
221  allocate(alpha(ip1jmp1,llm))
222  allocate(beta(ip1jmp1,llm))
223 
224  ! surface pressure
225  if (iflag_phys>2) then
226  ! specific value for CMIP5 aqua/terra planets
227  ! "Specify the initial dry mass to be equivalent to
228  ! a global mean surface pressure (101325 minus 245) Pa."
229  ps_glo(:)=101080.
230  else
231  ! use reference surface pressure
232  ps_glo(:)=preff
233  endif
234 
235  ! ground geopotential
236  phis_glo(:)=0.
237 
238  CALL pression( ip1jmp1, ap, bp, ps_glo, p )
239  if (pressure_exner) then
240  CALL exner_hyb( ip1jmp1, ps_glo, p,alpha,beta, pks, pk, pkf )
241  else
242  call exner_milieu(ip1jmp1,ps_glo,p,beta,pks,pk,pkf)
243  endif
244  CALL massdair(p,masse_glo)
245 
246  ! bulk initialization of temperature
247  teta_glo(:,:)=tetarappel(:,:)
248 
249  ! geopotential
250  CALL geopot(ip1jmp1,teta_glo,pk,pks,phis_glo,phi)
251 
252  ! winds
253  if (ok_geost) then
254  call ugeostr(phi,ucov_glo)
255  else
256  ucov_glo(:,:)=0.
257  endif
258  vcov(ijb_v:ije_v,1:llm)=0.
259 
260  ! bulk initialization of tracers
261  if (planet_type=="earth") then
262  ! Earth: first two tracers will be water
263  do i=1,nqtot
264  if (i == 1) q(ijb_u:ije_u,:,i)=1.e-10
265  if (i == 2) q(ijb_u:ije_u,:,i)=1.e-15
266  if (i.gt.2) q(ijb_u:ije_u,:,i)=0.
267  enddo
268  else
269  q(ijb_u:ije_u,:,:)=0
270  endif ! of if (planet_type=="earth")
271 
272  ! add random perturbation to temperature
273  idum = -1
274  zz = ran1(idum)
275  idum = 0
276  do l=1,llm
277  do ij=iip2,ip1jm
278  teta_glo(ij,l)=teta_glo(ij,l)*(1.+0.005*ran1(idum))
279  enddo
280  enddo
281 
282  ! maintain periodicity in longitude
283  do l=1,llm
284  do ij=1,ip1jmp1,iip1
285  teta_glo(ij+iim,l)=teta_glo(ij,l)
286  enddo
287  enddo
288 
289  ! copy data from global array to local array:
290  teta(ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:)
291  ucov(ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:)
292 ! vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:)
293  masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:)
294  ps(ijb_u:ije_u)=ps_glo(ijb_u:ije_u)
295  phis(ijb_u:ije_u)=phis_glo(ijb_u:ije_u)
296 
297  deallocate(teta_glo)
298  deallocate(ucov_glo)
299 ! deallocate(vcov_glo)
300  deallocate(masse_glo)
301  deallocate(ps_glo)
302  deallocate(phis_glo)
303  deallocate(alpha)
304  deallocate(beta)
305  ENDIF ! of IF (.NOT. read_start)
306  endif academic_case
307 
308 END SUBROUTINE iniacademic_loc