source: LMDZ5/branches/testing/libf/dyn3dmem/iniacademic_loc.F90 @ 2272

Last change on this file since 2272 was 2160, checked in by Laurent Fairhead, 10 years ago

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