source: LMDZ4/trunk/libf/dyn3dpar/iniacademic.F @ 1352

Last change on this file since 1352 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.1 KB
Line 
1!
2! $Header$
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
11c%W%    %G%
12c=======================================================================
13c
14c   Author:    Frederic Hourdin      original: 15/01/93
15c   -------
16c
17c   Subject:
18c   ------
19c
20c   Method:
21c   --------
22c
23c   Interface:
24c   ----------
25c
26c      Input:
27c      ------
28c
29c      Output:
30c      -------
31c
32c=======================================================================
33      IMPLICIT NONE
34c-----------------------------------------------------------------------
35c   Declararations:
36c   ---------------
37
38#include "dimensions.h"
39#include "paramet.h"
40#include "comvert.h"
41#include "comconst.h"
42#include "comgeom.h"
43#include "academic.h"
44#include "ener.h"
45#include "temps.h"
46#include "control.h"
47#include "iniprint.h"
48
49c   Arguments:
50c   ----------
51
52      real time_0
53
54c   variables dynamiques
55      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
56      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
57      REAL q(ip1jmp1,llm,nqtot)              ! champs advectes
58      REAL ps(ip1jmp1)                       ! pression  au sol
59      REAL masse(ip1jmp1,llm)                ! masse d'air
60      REAL phis(ip1jmp1)                     ! geopotentiel au sol
61
62c   Local:
63c   ------
64
65      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
66      REAL pks(ip1jmp1)                      ! exner au  sol
67      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
68      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
69      REAL phi(ip1jmp1,llm)                  ! geopotentiel
70      REAL ddsin,tetarappelj,tetarappell,zsig
71      real tetajl(jjp1,llm)
72      INTEGER i,j,l,lsup,ij
73
74      real zz,ran1
75      integer idum
76
77      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
78
79c-----------------------------------------------------------------------
80! 1. Initializations for Earth-like case
81! --------------------------------------
82      if (planet_type=="earth") then
83c
84        time_0=0.
85        day_ref=0
86        annee_ref=0
87
88        im         = iim
89        jm         = jjm
90        day_ini    = 0
91        omeg       = 4.*asin(1.)/86400.
92        rad    = 6371229.
93        g      = 9.8
94        daysec = 86400.
95        dtvr    = daysec/FLOAT(day_step)
96        zdtvr=dtvr
97        kappa  = 0.2857143
98        cpp    = 1004.70885
99        preff     = 101325.
100        pa        =  50000.
101        etot0      = 0.
102        ptot0      = 0.
103        ztot0      = 0.
104        stot0      = 0.
105        ang0       = 0.
106
107        CALL iniconst
108        CALL inigeom
109        CALL inifilr
110
111        ps=0.
112        phis=0.
113c---------------------------------------------------------------------
114
115        taurappel=10.*daysec
116
117c---------------------------------------------------------------------
118c   Calcul de la temperature potentielle :
119c   --------------------------------------
120
121        DO l=1,llm
122         zsig=ap(l)/preff+bp(l)
123         if (zsig.gt.0.3) then
124           lsup=l
125           tetarappell=1./8.*(-log(zsig)-.5)
126           DO j=1,jjp1
127             ddsin=sin(rlatu(j))-sin(pi/20.)
128             tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
129           ENDDO
130          else
131c   Choix isotherme au-dessus de 300 mbar
132           do j=1,jjp1
133             tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
134           enddo
135          endif ! of if (zsig.gt.0.3)
136        ENDDO ! of DO l=1,llm
137
138        do l=1,llm
139           do j=1,jjp1
140              do i=1,iip1
141                 ij=(j-1)*iip1+i
142                 tetarappel(ij,l)=tetajl(j,l)
143              enddo
144           enddo
145        enddo
146
147c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
148
149        ps=1.e5
150        phis=0.
151        CALL pression ( ip1jmp1, ap, bp, ps, p       )
152        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
153        CALL massdair(p,masse)
154
155c  intialisation du vent et de la temperature
156        teta(:,:)=tetarappel(:,:)
157        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
158        call ugeostr(phi,ucov)
159        vcov=0.
160        q(:,:,1   )=1.e-10
161        q(:,:,2   )=1.e-15
162        q(:,:,3:nqtot)=0.
163
164
165c   perturbation aleatoire sur la temperature
166        idum  = -1
167        zz = ran1(idum)
168        idum  = 0
169        do l=1,llm
170           do ij=iip2,ip1jm
171              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
172           enddo
173        enddo
174
175        do l=1,llm
176           do ij=1,ip1jmp1,iip1
177              teta(ij+iim,l)=teta(ij,l)
178           enddo
179        enddo
180
181
182
183c     PRINT *,' Appel test_period avec tetarappel '
184c     CALL  test_period ( ucov,vcov,tetarappel,q,p,phis )
185c     PRINT *,' Appel test_period avec teta '
186c     CALL  test_period ( ucov,vcov,teta,q,p,phis )
187
188c   initialisation d'un traceur sur une colonne
189        j=jjp1*3/4
190        i=iip1/2
191        ij=(j-1)*iip1+i
192        q(ij,:,3)=1.
193     
194      else
195        write(lunout,*)"iniacademic: planet types other than earth",
196     &                 " not implemented (yet)."
197        stop
198      endif ! of if (planet_type=="earth")
199      return
200      END
201c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.