source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3dpar/iniacademic.F90 @ 3693

Last change on this file since 3693 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.7 KB
Line 
1!
2! $Id: iniacademic.F90 1520 2011-05-23 11:37:09Z 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        ps(:)=preff
211        ! ground geopotential
212        phis(:)=0.
213
214        CALL pression ( ip1jmp1, ap, bp, ps, p       )
215        if (disvert_type.eq.1) then
216          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
217        elseif (disvert_type.eq.2) then
218          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
219        else
220          write(abort_message,*) "Wrong value for disvert_type: ", &
221                              disvert_type
222          call abort_gcm(modname,abort_message,0)
223        endif
224        CALL massdair(p,masse)
225
226        ! bulk initialization of temperature
227        teta(:,:)=tetarappel(:,:)
228
229        ! geopotential
230        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
231
232        ! winds
233        if (ok_geost) then
234           call ugeostr(phi,ucov)
235        else
236           ucov(:,:)=0.
237        endif
238        vcov(:,:)=0.
239
240        ! bulk initialization of tracers
241        if (planet_type=="earth") then
242           ! Earth: first two tracers will be water
243           do i=1,nqtot
244              if (i == 1) q(:,:,i)=1.e-10
245              if (i == 2) q(:,:,i)=1.e-15
246              if (i.gt.2) q(:,:,i)=0.
247           enddo
248        else
249           q(:,:,:)=0
250        endif ! of if (planet_type=="earth")
251
252        ! add random perturbation to temperature
253        idum  = -1
254        zz = ran1(idum)
255        idum  = 0
256        do l=1,llm
257           do ij=iip2,ip1jm
258              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
259           enddo
260        enddo
261
262        ! maintain periodicity in longitude
263        do l=1,llm
264           do ij=1,ip1jmp1,iip1
265              teta(ij+iim,l)=teta(ij,l)
266           enddo
267        enddo
268
269     ENDIF ! of IF (.NOT. read_start)
270  endif academic_case
271
272END SUBROUTINE iniacademic
Note: See TracBrowser for help on using the repository browser.