source: LMDZ5/branches/AI-cosp/libf/dyn3d/iniacademic.F90 @ 3738

Last change on this file since 3738 was 2270, checked in by crisi, 10 years ago

Adding isotopes in the dynamics and more generally tracers of tracers.
CRisi

  • 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: 9.1 KB
Line 
1!
2! $Id: iniacademic.F90 2270 2015-05-07 15:45:04Z adurocher $
3!
4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
5
6  USE filtreg_mod, ONLY: inifilr
7  USE infotrac
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
265              ! CRisi: init des isotopes
266              ! distill de Rayleigh très simplifiée
267              if (ok_isotopes) then
268                if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then         
269                   q(:,:,i)=q(:,:,iqpere(i))             &
270      &                  *tnat(iso_num(i))               &
271      &                  *(q(:,:,iqpere(i))/30.e-3)      &
272      &                  **(alpha_ideal(iso_num(i))-1)
273                endif               
274                if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then
275                  q(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i)))
276                endif
277              endif !if (ok_isotopes) then
278
279           enddo
280        else
281           q(:,:,:)=0
282        endif ! of if (planet_type=="earth")
283
284        if (ok_iso_verif) then
285           call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
286        endif !if (ok_iso_verif) then
287
288        ! add random perturbation to temperature
289        idum  = -1
290        zz = ran1(idum)
291        idum  = 0
292        do l=1,llm
293           do ij=iip2,ip1jm
294              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
295           enddo
296        enddo
297
298        ! maintain periodicity in longitude
299        do l=1,llm
300           do ij=1,ip1jmp1,iip1
301              teta(ij+iim,l)=teta(ij,l)
302           enddo
303        enddo
304
305     ENDIF ! of IF (.NOT. read_start)
306  endif academic_case
307
308END SUBROUTINE iniacademic
Note: See TracBrowser for help on using the repository browser.