source: LMDZ5/branches/LMDZ5-DOFOCO/libf/dyn3dmem/iniacademic_loc.F90 @ 4785

Last change on this file since 4785 was 1749, checked in by Ehouarn Millour, 12 years ago

Added handling of Newtonian (and Shallow Water) modes in dyn3dmem:

  • adapted calfis_loc.F and gcm.F
  • removed unused routines divgrad_p.F and gradiv_p.F
  • adapted iniacademic.F90 and sw_case_williamson.F to work in parallel and renamed them iniacademi_loc.F90 and sw_case_williamson_loc.F
  • fixed bug in exner_milieu_loc.F (filtreg_p form mod_filtreg_p should be used)

EM

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