source: LMDZ6/trunk/libf/dyn3d/iniacademic.F90 @ 5087

Last change on this file since 5087 was 5084, checked in by Laurent Fairhead, 12 months ago

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, 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
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.0 KB
Line 
1!
2! $Id: iniacademic.F90 5084 2024-07-19 16:40:44Z abarral $
3!
4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
5
6  USE filtreg_mod, ONLY: inifilr
7  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
8  USE control_mod, ONLY: day_step,planet_type
9  use exner_hyb_m, only: exner_hyb
10  use exner_milieu_m, only: exner_milieu
11#ifdef CPP_IOIPSL
12  USE IOIPSL, ONLY: getin
13#else
14  ! if not using IOIPSL, we still need to use (a local version of) getin
15  USE ioipsl_getincom, ONLY: getin
16#endif
17  USE Write_Field
18  USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm
19  USE logic_mod, ONLY: iflag_phys, read_start
20  USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
21  USE temps_mod, ONLY: annee_ref, day_ini, day_ref
22  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
23  USE readTracFiles_mod, ONLY: addPhase
24  use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
25  use netcdf, only : NF90_CLOSE, NF90_GET_VAR
26
27
28  !   Author:    Frederic Hourdin      original: 15/01/93
29  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
30  ! of the American Meteorological Society, 75, 1825.
31
32  IMPLICIT NONE
33
34  !   Declararations:
35  !   ---------------
36
37  include "dimensions.h"
38  include "paramet.h"
39  include "comgeom.h"
40  include "academic.h"
41  include "iniprint.h"
42
43  !   Arguments:
44  !   ----------
45
46  REAL,INTENT(OUT) :: time_0
47
48  !   fields
49  REAL,INTENT(OUT) :: vcov(ip1jm,llm) ! meridional covariant wind
50  REAL,INTENT(OUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
51  REAL,INTENT(OUT) :: teta(ip1jmp1,llm) ! potential temperature (K)
52  REAL,INTENT(OUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers (.../kg_of_air)
53  REAL,INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
54  REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass in grid cell (kg)
55  REAL,INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
56
57  !   Local:
58  !   ------
59
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, iq, iName, iPhase, iqParent
68
69  integer :: nid_relief,varid,ierr
70  real, dimension(iip1,jjp1) :: relief
71
72  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
73  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
74  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
75  LOGICAL ok_pv                ! Polar Vortex
76  REAL phi_pv,dphi_pv,gam_pv,tetanoise   ! Constantes pour polar vortex
77
78  real zz,ran1
79  integer idum
80
81  REAL zdtvr, tnat, alpha_ideal
82  LOGICAL,PARAMETER :: tnat1=.true.
83 
84  character(len=*),parameter :: modname="iniacademic"
85  character(len=80) :: abort_message
86
87  ! Sanity check: verify that options selected by user are not incompatible
88  if ((iflag_phys==1).and. .not. read_start) then
89    write(lunout,*) trim(modname)," error: if read_start is set to ", &
90    " false then iflag_phys should not be 1"
91    write(lunout,*) "You most likely want an aquaplanet initialisation", &
92    " (iflag_phys >= 100)"
93    call abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.false.",1)
94  endif
95 
96  !-----------------------------------------------------------------------
97  ! 1. Initializations for Earth-like case
98  ! --------------------------------------
99  !
100  ! initialize planet radius, rotation rate,...
101  call conf_planete
102
103  time_0=0.
104  day_ref=1
105  annee_ref=0
106
107  im         = iim
108  jm         = jjm
109  day_ini    = 1
110  dtvr    = daysec/REAL(day_step)
111  zdtvr=dtvr
112  etot0      = 0.
113  ptot0      = 0.
114  ztot0      = 0.
115  stot0      = 0.
116  ang0       = 0.
117
118  if (llm == 1) then
119     ! specific initializations for the shallow water case
120     kappa=1
121  endif
122
123  CALL iniconst
124  CALL inigeom
125  CALL inifilr
126
127
128  !------------------------------------------------------------------
129  ! Initialize pressure and mass field if read_start=.false.
130  !------------------------------------------------------------------
131
132  IF (.NOT. read_start) THEN
133
134     !------------------------------------------------------------------
135     ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
136     ! du modele.
137     ! On suppose que le fichier relief_in.nc est stoké sur une grille
138     ! iim*jjp1
139     ! Facile a créer à partir de la commande
140     ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
141     !------------------------------------------------------------------
142
143     relief=0.
144     ierr = NF90_OPEN ('relief_in.nc', NF90_NOWRITE,nid_relief)
145     if (ierr.EQ.NF90_NOERR) THEN
146         ierr=NF90_INQ_VARID(nid_relief,'RELIEF',varid)
147         if (ierr==NF90_NOERR) THEN
148              ierr=NF90_GET_VAR(nid_relief,varid,relief(1:iim,1:jjp1))
149              relief(iip1,:)=relief(1,:)
150         else
151              CALL abort_gcm ('iniacademic','variable RELIEF pas la',1)
152         endif
153     endif
154     ierr = NF90_CLOSE (nid_relief)
155
156     !------------------------------------------------------------------
157     ! Initialisation du geopotentiel au sol et de la pression
158     !------------------------------------------------------------------
159
160     print*,'relief=',minval(relief),maxval(relief),'g=',g
161     do j=1,jjp1
162        do i=1,iip1
163           phis((j-1)*iip1+i)=g*relief(i,j)
164        enddo
165     enddo
166     print*,'phis=',minval(phis),maxval(phis),'g=',g
167
168     ! ground geopotential
169     !phis(:)=0.
170     ps(:)=preff
171     CALL pression ( ip1jmp1, ap, bp, ps, p       )
172
173     if (pressure_exner) then
174       CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
175     else
176       call exner_milieu(ip1jmp1,ps,p,pks,pk)
177     endif
178     CALL massdair(p,masse)
179  ENDIF
180
181  if (llm == 1) then
182     ! initialize fields for the shallow water case, if required
183     if (.not.read_start) then
184        phis(:)=0.
185        q(:,:,:)=0
186        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
187     endif
188  endif
189
190  academic_case: if (iflag_phys >= 2) then
191     ! initializations
192
193     ! 1. local parameters
194     ! by convention, winter is in the southern hemisphere
195     ! Geostrophic wind or no wind?
196     ok_geost=.TRUE.
197     CALL getin('ok_geost',ok_geost)
198     ! Constants for Newtonian relaxation and friction
199     k_f=1.                !friction
200     CALL getin('k_j',k_f)
201     k_f=1./(daysec*k_f)
202     k_c_s=4.  !cooling surface
203     CALL getin('k_c_s',k_c_s)
204     k_c_s=1./(daysec*k_c_s)
205     k_c_a=40. !cooling free atm
206     CALL getin('k_c_a',k_c_a)
207     k_c_a=1./(daysec*k_c_a)
208     ! Constants for Teta equilibrium profile
209     teta0=315.     ! mean Teta (S.H. 315K)
210     CALL getin('teta0',teta0)
211     ttp=200.       ! Tropopause temperature (S.H. 200K)
212     CALL getin('ttp',ttp)
213     eps=0.         ! Deviation to N-S symmetry(~0-20K)
214     CALL getin('eps',eps)
215     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
216     CALL getin('delt_y',delt_y)
217     delt_z=10.     ! Vertical Gradient (S.H. 10K)
218     CALL getin('delt_z',delt_z)
219     ! Polar vortex
220     ok_pv=.false.
221     CALL getin('ok_pv',ok_pv)
222     phi_pv=-50.            ! Latitude of edge of vortex
223     CALL getin('phi_pv',phi_pv)
224     phi_pv=phi_pv*pi/180.
225     dphi_pv=5.             ! Width of the edge
226     CALL getin('dphi_pv',dphi_pv)
227     dphi_pv=dphi_pv*pi/180.
228     gam_pv=4.              ! -dT/dz vortex (in K/km)
229     CALL getin('gam_pv',gam_pv)
230     tetanoise=0.005
231     CALL getin('tetanoise',tetanoise)
232
233     ! 2. Initialize fields towards which to relax
234     ! Friction
235     knewt_g=k_c_a
236     DO l=1,llm
237        zsig=presnivs(l)/preff
238        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
239        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
240     ENDDO
241     DO j=1,jjp1
242        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
243     ENDDO
244
245     ! Potential temperature
246     DO l=1,llm
247        zsig=presnivs(l)/preff
248        tetastrat=ttp*zsig**(-kappa)
249        tetapv=tetastrat
250        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
251           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
252        ENDIF
253        DO j=1,jjp1
254           ! Troposphere
255           ddsin=sin(rlatu(j))
256           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
257                -delt_z*(1.-ddsin*ddsin)*log(zsig)
258           if (planet_type=="giant") then
259             tetajl(j,l)=teta0+(delt_y*                   &
260                ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
261                / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
262                -delt_z*log(zsig)
263           endif
264           ! Profil stratospherique isotherme (+vortex)
265           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
266           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
267           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
268        ENDDO
269     ENDDO
270
271     !          CALL writefield('theta_eq',tetajl)
272
273     do l=1,llm
274        do j=1,jjp1
275           do i=1,iip1
276              ij=(j-1)*iip1+i
277              tetarappel(ij,l)=tetajl(j,l)
278           enddo
279        enddo
280     enddo
281
282     ! 3. Initialize fields (if necessary)
283     IF (.NOT. read_start) THEN
284        ! bulk initialization of temperature
285        IF (iflag_phys>10000) THEN
286        ! Particular case to impose a constant temperature T0=0.01*iflag_physx
287           teta(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp)
288        ELSE
289           teta(:,:)=tetarappel(:,:)
290        ENDIF
291        ! geopotential
292        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
293
294        DO l=1,llm
295          print*,'presnivs,play,l',presnivs(l),(pk(1,l)/cpp)**(1./kappa)*preff
296         !pks(ij) = (cpp/preff) * ps(ij)
297         !pk(ij,1) = .5*pks(ij)
298         ! pk = cpp * (p/preff)^kappa
299        ENDDO
300
301        ! winds
302        if (ok_geost) then
303           call ugeostr(phi,ucov)
304        else
305           ucov(:,:)=0.
306        endif
307        vcov(:,:)=0.
308
309        ! bulk initialization of tracers
310        if (planet_type=="earth") then
311           ! Earth: first two tracers will be water
312           do iq=1,nqtot
313              q(:,:,iq)=0.
314              IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10
315              IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15
316
317              ! CRisi: init des isotopes
318              ! distill de Rayleigh très simplifiée
319              iName    = tracers(iq)%iso_iName
320              if (niso <= 0 .OR. iName <= 0) CYCLE
321              iPhase   = tracers(iq)%iso_iPhase
322              iqParent = tracers(iq)%iqParent
323              IF(tracers(iq)%iso_iZone == 0) THEN
324                 if (tnat1) then
325                         tnat=1.0
326                         alpha_ideal=1.0
327                         write(*,*) 'Attention dans iniacademic: alpha_ideal=1'
328                 else
329                    IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
330                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
331                 endif
332                 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
333              ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
334                 IF(tracers(iq)%iso_iZone == 1) THEN
335                    ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
336                    ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
337                    q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
338                 else !IF(tracers(iq)%iso_iZone == 1) THEN
339                    q(:,:,iq) = 0.
340                 endif !IF(tracers(iq)%iso_iZone == 1) THEN
341              END IF !IF(tracers(iq)%iso_iZone == 0) THEN
342           enddo
343        else
344           q(:,:,:)=0
345        endif ! of if (planet_type=="earth")
346
347        call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
348
349        ! add random perturbation to temperature
350        idum  = -1
351        zz = ran1(idum)
352        idum  = 0
353        do l=1,llm
354           do ij=iip2,ip1jm
355              teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum))
356           enddo
357        enddo
358
359        ! maintain periodicity in longitude
360        do l=1,llm
361           do ij=1,ip1jmp1,iip1
362              teta(ij+iim,l)=teta(ij,l)
363           enddo
364        enddo
365
366     ENDIF ! of IF (.NOT. read_start)
367  endif academic_case
368
369END SUBROUTINE iniacademic
Note: See TracBrowser for help on using the repository browser.