source: LMDZ5/branches/LMDZ5-DOFOCO/libf/dyn3d/iniacademic.F90 @ 3575

Last change on this file since 3575 was 1625, checked in by lguez, 12 years ago

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