source: LMDZ6/branches/SETHET_DECOUPLE/libf/dyn3dmem/iniacademic_loc.F90 @ 5437

Last change on this file since 5437 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

  • as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

  • due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

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