source: LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/iniacademic.F @ 5405

Last change on this file since 5405 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.4 KB
Line 
1!
2! $Id: iniacademic.F 1446 2010-10-22 09:27:25Z abarral $
3!
4c
5c
6      SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
7
8      USE filtreg_mod
9      USE infotrac, ONLY : nqtot
10      USE control_mod, ONLY: day_step,planet_type
11#ifdef CPP_IOIPSL
12      USE IOIPSL
13#else
14! if not using IOIPSL, we still need to use (a local version of) getin
15      USE ioipsl_getincom
16#endif
17      USE Write_Field
18 
19
20c%W%    %G%
21c=======================================================================
22c
23c   Author:    Frederic Hourdin      original: 15/01/93
24c   -------
25c
26c   Subject:
27c   ------
28c
29c   Method:
30c   --------
31c
32c   Interface:
33c   ----------
34c
35c      Input:
36c      ------
37c
38c      Output:
39c      -------
40c
41c=======================================================================
42      IMPLICIT NONE
43c-----------------------------------------------------------------------
44c   Declararations:
45c   ---------------
46
47#include "dimensions.h"
48#include "paramet.h"
49#include "comvert.h"
50#include "comconst.h"
51#include "comgeom.h"
52#include "academic.h"
53#include "ener.h"
54#include "temps.h"
55#include "iniprint.h"
56#include "logic.h"
57
58c   Arguments:
59c   ----------
60
61      real time_0
62
63c   variables dynamiques
64      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
65      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
66      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
67      REAL ps(ip1jmp1)                       ! pression  au sol
68      REAL masse(ip1jmp1,llm)                ! masse d'air
69      REAL phis(ip1jmp1)                     ! geopotentiel au sol
70
71c   Local:
72c   ------
73
74      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
75      REAL pks(ip1jmp1)                      ! exner au  sol
76      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
77      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
78      REAL phi(ip1jmp1,llm)                  ! geopotentiel
79      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
80      real tetajl(jjp1,llm)
81      INTEGER i,j,l,lsup,ij
82
83      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
84      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
85      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
86      LOGICAL ok_pv                ! Polar Vortex
87      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
88     
89      real zz,ran1
90      integer idum
91
92      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
93
94c-----------------------------------------------------------------------
95! 1. Initializations for Earth-like case
96! --------------------------------------
97c
98        ! initialize planet radius, rotation rate,...
99        call conf_planete
100
101        time_0=0.
102        day_ref=1
103        annee_ref=0
104
105        im         = iim
106        jm         = jjm
107        day_ini    = 1
108        dtvr    = daysec/REAL(day_step)
109        zdtvr=dtvr
110        etot0      = 0.
111        ptot0      = 0.
112        ztot0      = 0.
113        stot0      = 0.
114        ang0       = 0.
115
116        if (llm.eq.1) then
117          ! specific initializations for the shallow water case
118          kappa=1
119        endif
120       
121        CALL iniconst
122        CALL inigeom
123        CALL inifilr
124
125        if (llm.eq.1) then
126          ! initialize fields for the shallow water case, if required
127          if (.not.read_start) then
128            phis(:)=0.
129            q(:,:,:)=0
130            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
131          endif
132        endif
133
134        if (iflag_phys.eq.2) then
135          ! initializations for the academic case
136         
137!         if (planet_type=="earth") then
138
139          ! 1. local parameters
140          ! by convention, winter is in the southern hemisphere
141          ! Geostrophic wind or no wind?
142          ok_geost=.TRUE.
143          CALL getin('ok_geost',ok_geost)
144          ! Constants for Newtonian relaxation and friction
145          k_f=1.                !friction
146          CALL getin('k_j',k_f)
147          k_f=1./(daysec*k_f)
148          k_c_s=4.  !cooling surface
149          CALL getin('k_c_s',k_c_s)
150          k_c_s=1./(daysec*k_c_s)
151          k_c_a=40. !cooling free atm
152          CALL getin('k_c_a',k_c_a)
153          k_c_a=1./(daysec*k_c_a)
154          ! Constants for Teta equilibrium profile
155          teta0=315.     ! mean Teta (S.H. 315K)
156          CALL getin('teta0',teta0)
157          ttp=200.       ! Tropopause temperature (S.H. 200K)
158          CALL getin('ttp',ttp)
159          eps=0.         ! Deviation to N-S symmetry(~0-20K)
160          CALL getin('eps',eps)
161          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
162          CALL getin('delt_y',delt_y)
163          delt_z=10.     ! Vertical Gradient (S.H. 10K)
164          CALL getin('delt_z',delt_z)
165          ! Polar vortex
166          ok_pv=.false.
167          CALL getin('ok_pv',ok_pv)
168          phi_pv=-50.            ! Latitude of edge of vortex
169          CALL getin('phi_pv',phi_pv)
170          phi_pv=phi_pv*pi/180.
171          dphi_pv=5.             ! Width of the edge
172          CALL getin('dphi_pv',dphi_pv)
173          dphi_pv=dphi_pv*pi/180.
174          gam_pv=4.              ! -dT/dz vortex (in K/km)
175          CALL getin('gam_pv',gam_pv)
176         
177          ! 2. Initialize fields towards which to relax
178          ! Friction
179          knewt_g=k_c_a
180          DO l=1,llm
181           zsig=presnivs(l)/preff
182           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
183           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
184          ENDDO
185          DO j=1,jjp1
186            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
187          ENDDO
188         
189          ! Potential temperature
190          DO l=1,llm
191           zsig=presnivs(l)/preff
192           tetastrat=ttp*zsig**(-kappa)
193           tetapv=tetastrat
194           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
195               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
196           ENDIF
197           DO j=1,jjp1
198             ! Troposphere
199             ddsin=sin(rlatu(j))
200             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
201     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
202             ! Profil stratospherique isotherme (+vortex)
203             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
204             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
205             tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
206           ENDDO
207          ENDDO ! of DO l=1,llm
208 
209!          CALL writefield('theta_eq',tetajl)
210
211          do l=1,llm
212            do j=1,jjp1
213              do i=1,iip1
214                 ij=(j-1)*iip1+i
215                 tetarappel(ij,l)=tetajl(j,l)
216              enddo
217            enddo
218          enddo
219
220
221!         else
222!          write(lunout,*)"iniacademic: planet types other than earth",
223!     &                   " not implemented (yet)."
224!          stop
225!         endif ! of if (planet_type=="earth")
226
227          ! 3. Initialize fields (if necessary)
228          IF (.NOT. read_start) THEN
229            ! surface pressure
230            ps(:)=preff
231            ! ground geopotential
232            phis(:)=0.
233           
234            CALL pression ( ip1jmp1, ap, bp, ps, p       )
235            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
236            CALL massdair(p,masse)
237
238            ! bulk initialization of temperature
239            teta(:,:)=tetarappel(:,:)
240           
241            ! geopotential
242            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
243           
244            ! winds
245            if (ok_geost) then
246              call ugeostr(phi,ucov)
247            else
248              ucov(:,:)=0.
249            endif
250            vcov(:,:)=0.
251           
252            ! bulk initialization of tracers
253            if (planet_type=="earth") then
254              ! Earth: first two tracers will be water
255              do i=1,nqtot
256                if (i.eq.1) q(:,:,i)=1.e-10
257                if (i.eq.2) q(:,:,i)=1.e-15
258                if (i.gt.2) q(:,:,i)=0.
259              enddo
260            else
261              q(:,:,:)=0
262            endif ! of if (planet_type=="earth")
263
264            ! add random perturbation to temperature
265            idum  = -1
266            zz = ran1(idum)
267            idum  = 0
268            do l=1,llm
269              do ij=iip2,ip1jm
270                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
271              enddo
272            enddo
273
274            ! maintain periodicity in longitude
275            do l=1,llm
276              do ij=1,ip1jmp1,iip1
277                teta(ij+iim,l)=teta(ij,l)
278              enddo
279            enddo
280
281          ENDIF ! of IF (.NOT. read_start)
282        endif ! of if (iflag_phys.eq.2)
283       
284      END
285c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.