source: LMDZ5/branches/testing/libf/dyn3d/iniacademic.F90 @ 2056

Last change on this file since 2056 was 2056, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes r1997:2055 into testing branch

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