source: trunk/libf/dyn3dpar/iniacademic.F @ 7

Last change on this file since 7 was 7, checked in by emillour, 14 years ago

Mise a niveau de la dynamique par rapport a la version 1447 de LMDZ5
Voir les details dans chantiers/commit_v7.log
File size: 8.4 KB
RevLine 
[1]1!
[7]2! $Id: iniacademic.F 1446 2010-10-22 09:27:25Z emillour $
[1]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
[7]10      USE control_mod, ONLY: day_step,planet_type
[1]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.
[7]129            q(:,:,:)=0
[1]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         
[7]137!         if (planet_type=="earth") then
138
[1]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
[7]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
[1]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
[7]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")
[1]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
[7]274            ! maintain periodicity in longitude
[1]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.