source: LMDZ5/trunk/libf/dyn3d/iniacademic.F90 @ 2119

Last change on this file since 2119 was 2087, checked in by fhourdin, 10 years ago

Correction de bug pour gfortran qui refuse if (rerad_start==.false.)
change en if (.not. read_start)
  • 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: 8.3 KB
Line 
1!
2! $Id: iniacademic.F90 2087 2014-07-10 07:18:06Z fhourdin $
3!
4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
5
6  USE filtreg_mod, ONLY: inifilr
7  USE infotrac, ONLY : nqtot
8  USE control_mod, ONLY: day_step,planet_type
9#ifdef CPP_IOIPSL
10  USE IOIPSL, ONLY: getin
11#else
12  ! if not using IOIPSL, we still need to use (a local version of) getin
13  USE ioipsl_getincom, ONLY: getin
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,INTENT(OUT) :: time_0
43
44  !   fields
45  REAL,INTENT(OUT) :: vcov(ip1jm,llm) ! meridional covariant wind
46  REAL,INTENT(OUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
47  REAL,INTENT(OUT) :: teta(ip1jmp1,llm) ! potential temperature (K)
48  REAL,INTENT(OUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers (.../kg_of_air)
49  REAL,INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
50  REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass in grid cell (kg)
51  REAL,INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
52
53  !   Local:
54  !   ------
55
56  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
57  REAL pks(ip1jmp1)                      ! exner au  sol
58  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
59  REAL phi(ip1jmp1,llm)                  ! geopotentiel
60  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
61  real tetastrat ! potential temperature in the stratosphere, in K
62  real tetajl(jjp1,llm)
63  INTEGER i,j,l,lsup,ij
64
65  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
66  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
67  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
68  LOGICAL ok_pv                ! Polar Vortex
69  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
70
71  real zz,ran1
72  integer idum
73
74  REAL zdtvr
75 
76  character(len=*),parameter :: modname="iniacademic"
77  character(len=80) :: abort_message
78
79
80  ! Sanity check: verify that options selected by user are not incompatible
81  if ((iflag_phys==1).and. .not. read_start) then
82    write(lunout,*) trim(modname)," error: if read_start is set to ", &
83    " false then iflag_phys should not be 1"
84    write(lunout,*) "You most likely want an aquaplanet initialisation", &
85    " (iflag_phys >= 100)"
86    call abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.false.",1)
87  endif
88 
89  !-----------------------------------------------------------------------
90  ! 1. Initializations for Earth-like case
91  ! --------------------------------------
92  !
93  ! initialize planet radius, rotation rate,...
94  call conf_planete
95
96  time_0=0.
97  day_ref=1
98  annee_ref=0
99
100  im         = iim
101  jm         = jjm
102  day_ini    = 1
103  dtvr    = daysec/REAL(day_step)
104  zdtvr=dtvr
105  etot0      = 0.
106  ptot0      = 0.
107  ztot0      = 0.
108  stot0      = 0.
109  ang0       = 0.
110
111  if (llm == 1) then
112     ! specific initializations for the shallow water case
113     kappa=1
114  endif
115
116  CALL iniconst
117  CALL inigeom
118  CALL inifilr
119
120  if (llm == 1) then
121     ! initialize fields for the shallow water case, if required
122     if (.not.read_start) then
123        phis(:)=0.
124        q(:,:,:)=0
125        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
126     endif
127  endif
128
129  academic_case: if (iflag_phys >= 2) then
130     ! initializations
131
132     ! 1. local parameters
133     ! by convention, winter is in the southern hemisphere
134     ! Geostrophic wind or no wind?
135     ok_geost=.TRUE.
136     CALL getin('ok_geost',ok_geost)
137     ! Constants for Newtonian relaxation and friction
138     k_f=1.                !friction
139     CALL getin('k_j',k_f)
140     k_f=1./(daysec*k_f)
141     k_c_s=4.  !cooling surface
142     CALL getin('k_c_s',k_c_s)
143     k_c_s=1./(daysec*k_c_s)
144     k_c_a=40. !cooling free atm
145     CALL getin('k_c_a',k_c_a)
146     k_c_a=1./(daysec*k_c_a)
147     ! Constants for Teta equilibrium profile
148     teta0=315.     ! mean Teta (S.H. 315K)
149     CALL getin('teta0',teta0)
150     ttp=200.       ! Tropopause temperature (S.H. 200K)
151     CALL getin('ttp',ttp)
152     eps=0.         ! Deviation to N-S symmetry(~0-20K)
153     CALL getin('eps',eps)
154     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
155     CALL getin('delt_y',delt_y)
156     delt_z=10.     ! Vertical Gradient (S.H. 10K)
157     CALL getin('delt_z',delt_z)
158     ! Polar vortex
159     ok_pv=.false.
160     CALL getin('ok_pv',ok_pv)
161     phi_pv=-50.            ! Latitude of edge of vortex
162     CALL getin('phi_pv',phi_pv)
163     phi_pv=phi_pv*pi/180.
164     dphi_pv=5.             ! Width of the edge
165     CALL getin('dphi_pv',dphi_pv)
166     dphi_pv=dphi_pv*pi/180.
167     gam_pv=4.              ! -dT/dz vortex (in K/km)
168     CALL getin('gam_pv',gam_pv)
169
170     ! 2. Initialize fields towards which to relax
171     ! Friction
172     knewt_g=k_c_a
173     DO l=1,llm
174        zsig=presnivs(l)/preff
175        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
176        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
177     ENDDO
178     DO j=1,jjp1
179        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
180     ENDDO
181
182     ! Potential temperature
183     DO l=1,llm
184        zsig=presnivs(l)/preff
185        tetastrat=ttp*zsig**(-kappa)
186        tetapv=tetastrat
187        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
188           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
189        ENDIF
190        DO j=1,jjp1
191           ! Troposphere
192           ddsin=sin(rlatu(j))
193           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
194                -delt_z*(1.-ddsin*ddsin)*log(zsig)
195           if (planet_type=="giant") then
196             tetajl(j,l)=teta0+(delt_y*                   &
197                ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
198                / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
199                -delt_z*log(zsig)
200           endif
201           ! Profil stratospherique isotherme (+vortex)
202           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
203           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
204           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
205        ENDDO
206     ENDDO
207
208     !          CALL writefield('theta_eq',tetajl)
209
210     do l=1,llm
211        do j=1,jjp1
212           do i=1,iip1
213              ij=(j-1)*iip1+i
214              tetarappel(ij,l)=tetajl(j,l)
215           enddo
216        enddo
217     enddo
218
219     ! 3. Initialize fields (if necessary)
220     IF (.NOT. read_start) THEN
221        ! surface pressure
222        if (iflag_phys>2) then
223           ! specific value for CMIP5 aqua/terra planets
224           ! "Specify the initial dry mass to be equivalent to
225           !  a global mean surface pressure (101325 minus 245) Pa."
226           ps(:)=101080. 
227        else
228           ! use reference surface pressure
229           ps(:)=preff
230        endif
231       
232        ! ground geopotential
233        phis(:)=0.
234
235        CALL pression ( ip1jmp1, ap, bp, ps, p       )
236        if (pressure_exner) then
237          CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
238        else
239          call exner_milieu(ip1jmp1,ps,p,pks,pk)
240        endif
241        CALL massdair(p,masse)
242
243        ! bulk initialization of temperature
244        teta(:,:)=tetarappel(:,:)
245
246        ! geopotential
247        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
248
249        ! winds
250        if (ok_geost) then
251           call ugeostr(phi,ucov)
252        else
253           ucov(:,:)=0.
254        endif
255        vcov(:,:)=0.
256
257        ! bulk initialization of tracers
258        if (planet_type=="earth") then
259           ! Earth: first two tracers will be water
260           do i=1,nqtot
261              if (i == 1) q(:,:,i)=1.e-10
262              if (i == 2) q(:,:,i)=1.e-15
263              if (i.gt.2) q(:,:,i)=0.
264           enddo
265        else
266           q(:,:,:)=0
267        endif ! of if (planet_type=="earth")
268
269        ! add random perturbation to temperature
270        idum  = -1
271        zz = ran1(idum)
272        idum  = 0
273        do l=1,llm
274           do ij=iip2,ip1jm
275              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
276           enddo
277        enddo
278
279        ! maintain periodicity in longitude
280        do l=1,llm
281           do ij=1,ip1jmp1,iip1
282              teta(ij+iim,l)=teta(ij,l)
283           enddo
284        enddo
285
286     ENDIF ! of IF (.NOT. read_start)
287  endif academic_case
288
289END SUBROUTINE iniacademic
Note: See TracBrowser for help on using the repository browser.