source: trunk/libf/dyn3d/iniacademic.F90 @ 125

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

EM: suite mise au point discretisation verticale et quelques corrections de bugs dans le version de reference parallele.

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