source: LMDZ5/branches/LMDZ5-dev032011/libf/dyn3dpar/iniacademic.F90 @ 4411

Last change on this file since 4411 was 1492, checked in by Laurent Fairhead, 14 years ago

Merge of development branch LMDZ5V2.0-dev r1455:r1491 into the trunk.
Validation made locally: restart files are strictly equal between the HEAD of the trunk
and r1491 of LMDZ5V2.0-dev


Synchro de la branche de développement LMDZ5V2.0-dev r1455:r1491 et de la trunk
Validation faite en local: les fichiers restart sont équivalents entre la HEAD de la trunk
et la révision r1491 de LMDZ5V2.0-dev

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.0 KB
Line 
1!
2! $Id: iniacademic.F90 1492 2011-03-08 08:10:25Z lguez $
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  !-----------------------------------------------------------------------
75  ! 1. Initializations for Earth-like case
76  ! --------------------------------------
77  !
78  ! initialize planet radius, rotation rate,...
79  call conf_planete
80
81  time_0=0.
82  day_ref=1
83  annee_ref=0
84
85  im         = iim
86  jm         = jjm
87  day_ini    = 1
88  dtvr    = daysec/REAL(day_step)
89  zdtvr=dtvr
90  etot0      = 0.
91  ptot0      = 0.
92  ztot0      = 0.
93  stot0      = 0.
94  ang0       = 0.
95
96  if (llm == 1) then
97     ! specific initializations for the shallow water case
98     kappa=1
99  endif
100
101  CALL iniconst
102  CALL inigeom
103  CALL inifilr
104
105  if (llm == 1) then
106     ! initialize fields for the shallow water case, if required
107     if (.not.read_start) then
108        phis(:)=0.
109        q(:,:,:)=0
110        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
111     endif
112  endif
113
114  academic_case: if (iflag_phys == 2) then
115     ! initializations
116
117     ! 1. local parameters
118     ! by convention, winter is in the southern hemisphere
119     ! Geostrophic wind or no wind?
120     ok_geost=.TRUE.
121     CALL getin('ok_geost',ok_geost)
122     ! Constants for Newtonian relaxation and friction
123     k_f=1.                !friction
124     CALL getin('k_j',k_f)
125     k_f=1./(daysec*k_f)
126     k_c_s=4.  !cooling surface
127     CALL getin('k_c_s',k_c_s)
128     k_c_s=1./(daysec*k_c_s)
129     k_c_a=40. !cooling free atm
130     CALL getin('k_c_a',k_c_a)
131     k_c_a=1./(daysec*k_c_a)
132     ! Constants for Teta equilibrium profile
133     teta0=315.     ! mean Teta (S.H. 315K)
134     CALL getin('teta0',teta0)
135     ttp=200.       ! Tropopause temperature (S.H. 200K)
136     CALL getin('ttp',ttp)
137     eps=0.         ! Deviation to N-S symmetry(~0-20K)
138     CALL getin('eps',eps)
139     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
140     CALL getin('delt_y',delt_y)
141     delt_z=10.     ! Vertical Gradient (S.H. 10K)
142     CALL getin('delt_z',delt_z)
143     ! Polar vortex
144     ok_pv=.false.
145     CALL getin('ok_pv',ok_pv)
146     phi_pv=-50.            ! Latitude of edge of vortex
147     CALL getin('phi_pv',phi_pv)
148     phi_pv=phi_pv*pi/180.
149     dphi_pv=5.             ! Width of the edge
150     CALL getin('dphi_pv',dphi_pv)
151     dphi_pv=dphi_pv*pi/180.
152     gam_pv=4.              ! -dT/dz vortex (in K/km)
153     CALL getin('gam_pv',gam_pv)
154
155     ! 2. Initialize fields towards which to relax
156     ! Friction
157     knewt_g=k_c_a
158     DO l=1,llm
159        zsig=presnivs(l)/preff
160        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
161        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
162     ENDDO
163     DO j=1,jjp1
164        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
165     ENDDO
166
167     ! Potential temperature
168     DO l=1,llm
169        zsig=presnivs(l)/preff
170        tetastrat=ttp*zsig**(-kappa)
171        tetapv=tetastrat
172        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
173           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
174        ENDIF
175        DO j=1,jjp1
176           ! Troposphere
177           ddsin=sin(rlatu(j))
178           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
179                -delt_z*(1.-ddsin*ddsin)*log(zsig)
180           ! Profil stratospherique isotherme (+vortex)
181           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
182           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
183           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
184        ENDDO
185     ENDDO
186
187     !          CALL writefield('theta_eq',tetajl)
188
189     do l=1,llm
190        do j=1,jjp1
191           do i=1,iip1
192              ij=(j-1)*iip1+i
193              tetarappel(ij,l)=tetajl(j,l)
194           enddo
195        enddo
196     enddo
197
198     ! 3. Initialize fields (if necessary)
199     IF (.NOT. read_start) THEN
200        ! surface pressure
201        ps(:)=preff
202        ! ground geopotential
203        phis(:)=0.
204
205        CALL pression ( ip1jmp1, ap, bp, ps, p       )
206        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
207        CALL massdair(p,masse)
208
209        ! bulk initialization of temperature
210        teta(:,:)=tetarappel(:,:)
211
212        ! geopotential
213        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
214
215        ! winds
216        if (ok_geost) then
217           call ugeostr(phi,ucov)
218        else
219           ucov(:,:)=0.
220        endif
221        vcov(:,:)=0.
222
223        ! bulk initialization of tracers
224        if (planet_type=="earth") then
225           ! Earth: first two tracers will be water
226           do i=1,nqtot
227              if (i == 1) q(:,:,i)=1.e-10
228              if (i == 2) q(:,:,i)=1.e-15
229              if (i.gt.2) q(:,:,i)=0.
230           enddo
231        else
232           q(:,:,:)=0
233        endif ! of if (planet_type=="earth")
234
235        ! add random perturbation to temperature
236        idum  = -1
237        zz = ran1(idum)
238        idum  = 0
239        do l=1,llm
240           do ij=iip2,ip1jm
241              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
242           enddo
243        enddo
244
245        ! maintain periodicity in longitude
246        do l=1,llm
247           do ij=1,ip1jmp1,iip1
248              teta(ij+iim,l)=teta(ij,l)
249           enddo
250        enddo
251
252     ENDIF ! of IF (.NOT. read_start)
253  endif academic_case
254
255END SUBROUTINE iniacademic
Note: See TracBrowser for help on using the repository browser.