source: LMDZ6/branches/contrails/libf/dyn3d/iniacademic.f90 @ 5440

Last change on this file since 5440 was 5292, checked in by abarral, 8 weeks ago

Move academic.h chem.h chem_spla.h to module
  • 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.4 KB
Line 
1!
2! $Id: iniacademic.f90 5292 2024-10-28 15:58:32Z evignon $
3!
4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
5
6USE academic_mod_h
7    USE iniprint_mod_h
8  USE comgeom_mod_h
9  USE filtreg_mod, ONLY: inifilr
10  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase
11  USE control_mod, ONLY: day_step,planet_type
12  use exner_hyb_m, only: exner_hyb
13  use exner_milieu_m, only: exner_milieu
14  USE IOIPSL, ONLY: getin
15
16  USE Write_Field
17  USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm
18  USE logic_mod, ONLY: iflag_phys, read_start
19  USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
20  USE temps_mod, ONLY: annee_ref, day_ini, day_ref
21  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
22  use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
23  use netcdf, only : NF90_CLOSE, NF90_GET_VAR
24  USE iso_params_mod   ! tnat_* and alpha_ideal_*
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  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
32USE paramet_mod_h
33IMPLICIT NONE
34
35  !   Declararations:
36  !   ---------------
37
38
39
40
41  !   Arguments:
42  !   ----------
43
44  REAL,INTENT(OUT) :: time_0
45
46  !   fields
47  REAL,INTENT(OUT) :: vcov(ip1jm,llm) ! meridional covariant wind
48  REAL,INTENT(OUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
49  REAL,INTENT(OUT) :: teta(ip1jmp1,llm) ! potential temperature (K)
50  REAL,INTENT(OUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers (.../kg_of_air)
51  REAL,INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
52  REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass in grid cell (kg)
53  REAL,INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
54
55  !   Local:
56  !   ------
57
58  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
59  REAL pks(ip1jmp1)                      ! exner au  sol
60  REAL pk(ip1jmp1,llm)                   ! exner 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, iq, iName, iPhase, iqParent
66
67  integer :: nid_relief,varid,ierr
68  real, dimension(iip1,jjp1) :: relief
69
70  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
71  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
72  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
73  LOGICAL ok_pv                ! Polar Vortex
74  REAL phi_pv,dphi_pv,gam_pv,tetanoise   ! Constantes pour polar vortex
75
76  real zz,ran1
77  integer idum
78
79  REAL zdtvr, tnat, alpha_ideal
80  LOGICAL :: ltnat1
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
126  !------------------------------------------------------------------
127  ! Initialize pressure and mass field if read_start=.false.
128  !------------------------------------------------------------------
129
130  IF (.NOT. read_start) THEN
131
132     !------------------------------------------------------------------
133     ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
134     ! du modele.
135     ! On suppose que le fichier relief_in.nc est stoké sur une grille
136     ! iim*jjp1
137     ! Facile a créer à partir de la commande
138     ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
139     !------------------------------------------------------------------
140
141     relief=0.
142     ierr = NF90_OPEN ('relief_in.nc', NF90_NOWRITE,nid_relief)
143     if (ierr.EQ.NF90_NOERR) THEN
144         ierr=NF90_INQ_VARID(nid_relief,'RELIEF',varid)
145         if (ierr==NF90_NOERR) THEN
146              ierr=NF90_GET_VAR(nid_relief,varid,relief(1:iim,1:jjp1))
147              relief(iip1,:)=relief(1,:)
148         else
149              CALL abort_gcm ('iniacademic','variable RELIEF pas la',1)
150         endif
151     endif
152     ierr = NF90_CLOSE (nid_relief)
153
154     !------------------------------------------------------------------
155     ! Initialisation du geopotentiel au sol et de la pression
156     !------------------------------------------------------------------
157
158     print*,'relief=',minval(relief),maxval(relief),'g=',g
159     do j=1,jjp1
160        do i=1,iip1
161           phis((j-1)*iip1+i)=g*relief(i,j)
162        enddo
163     enddo
164     print*,'phis=',minval(phis),maxval(phis),'g=',g
165
166     ! ground geopotential
167     !phis(:)=0.
168     ps(:)=preff
169     CALL pression ( ip1jmp1, ap, bp, ps, p       )
170
171     if (pressure_exner) then
172       CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
173     else
174       call exner_milieu(ip1jmp1,ps,p,pks,pk)
175     endif
176     CALL massdair(p,masse)
177  ENDIF
178
179  if (llm == 1) then
180     ! initialize fields for the shallow water case, if required
181     if (.not.read_start) then
182        phis(:)=0.
183        q(:,:,:)=0
184        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
185     endif
186  endif
187
188  academic_case: if (iflag_phys >= 2) then
189     ! initializations
190
191     ! 1. local parameters
192     ! by convention, winter is in the southern hemisphere
193     ! Geostrophic wind or no wind?
194     ok_geost=.TRUE.
195     CALL getin('ok_geost',ok_geost)
196     ! Constants for Newtonian relaxation and friction
197     k_f=1.                !friction
198     CALL getin('k_j',k_f)
199     k_f=1./(daysec*k_f)
200     k_c_s=4.  !cooling surface
201     CALL getin('k_c_s',k_c_s)
202     k_c_s=1./(daysec*k_c_s)
203     k_c_a=40. !cooling free atm
204     CALL getin('k_c_a',k_c_a)
205     k_c_a=1./(daysec*k_c_a)
206     ! Constants for Teta equilibrium profile
207     teta0=315.     ! mean Teta (S.H. 315K)
208     CALL getin('teta0',teta0)
209     ttp=200.       ! Tropopause temperature (S.H. 200K)
210     CALL getin('ttp',ttp)
211     eps=0.         ! Deviation to N-S symmetry(~0-20K)
212     CALL getin('eps',eps)
213     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
214     CALL getin('delt_y',delt_y)
215     delt_z=10.     ! Vertical Gradient (S.H. 10K)
216     CALL getin('delt_z',delt_z)
217     ! Polar vortex
218     ok_pv=.false.
219     CALL getin('ok_pv',ok_pv)
220     phi_pv=-50.            ! Latitude of edge of vortex
221     CALL getin('phi_pv',phi_pv)
222     phi_pv=phi_pv*pi/180.
223     dphi_pv=5.             ! Width of the edge
224     CALL getin('dphi_pv',dphi_pv)
225     dphi_pv=dphi_pv*pi/180.
226     gam_pv=4.              ! -dT/dz vortex (in K/km)
227     CALL getin('gam_pv',gam_pv)
228     tetanoise=0.005
229     CALL getin('tetanoise',tetanoise)
230
231     ! 2. Initialize fields towards which to relax
232     ! Friction
233     knewt_g=k_c_a
234     DO l=1,llm
235        zsig=presnivs(l)/preff
236        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
237        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
238     ENDDO
239     DO j=1,jjp1
240        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
241     ENDDO
242
243     ! Potential temperature
244     DO l=1,llm
245        zsig=presnivs(l)/preff
246        tetastrat=ttp*zsig**(-kappa)
247        tetapv=tetastrat
248        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
249           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
250        ENDIF
251        DO j=1,jjp1
252           ! Troposphere
253           ddsin=sin(rlatu(j))
254           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
255                -delt_z*(1.-ddsin*ddsin)*log(zsig)
256           if (planet_type=="giant") then
257             tetajl(j,l)=teta0+(delt_y*                   &
258                ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
259                / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
260                -delt_z*log(zsig)
261           endif
262           ! Profil stratospherique isotherme (+vortex)
263           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
264           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
265           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
266        ENDDO
267     ENDDO
268
269     !          CALL writefield('theta_eq',tetajl)
270
271     do l=1,llm
272        do j=1,jjp1
273           do i=1,iip1
274              ij=(j-1)*iip1+i
275              tetarappel(ij,l)=tetajl(j,l)
276           enddo
277        enddo
278     enddo
279
280     ! 3. Initialize fields (if necessary)
281     IF (.NOT. read_start) THEN
282        ! bulk initialization of temperature
283        IF (iflag_phys>10000) THEN
284        ! Particular case to impose a constant temperature T0=0.01*iflag_physx
285           teta(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp)
286        ELSE
287           teta(:,:)=tetarappel(:,:)
288        ENDIF
289        ! geopotential
290        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
291
292        DO l=1,llm
293          print*,'presnivs,play,l',presnivs(l),(pk(1,l)/cpp)**(1./kappa)*preff
294         !pks(ij) = (cpp/preff) * ps(ij)
295         !pk(ij,1) = .5*pks(ij)
296         ! pk = cpp * (p/preff)^kappa
297        ENDDO
298
299        ! winds
300        if (ok_geost) then
301           call ugeostr(phi,ucov)
302        else
303           ucov(:,:)=0.
304        endif
305        vcov(:,:)=0.
306
307        ! bulk initialization of tracers
308        if (planet_type=="earth") then
309           ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
310           ! Earth: first two tracers will be water
311           do iq=1,nqtot
312              q(:,:,iq)=0.
313              IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10
314              IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15
315
316              ! CRisi: init des isotopes
317              ! distill de Rayleigh très simplifiée
318              iName    = tracers(iq)%iso_iName
319              if (niso <= 0 .OR. iName <= 0) CYCLE
320              iPhase   = tracers(iq)%iso_iPhase
321              iqParent = tracers(iq)%iqParent
322              IF(tracers(iq)%iso_iZone == 0) THEN
323                 IF(ltnat1) THEN
324                    tnat = 1.0
325                    alpha_ideal = 1.0
326                    WRITE(lunout, *)'In '//TRIM(modname)//': !!!  Beware: alpha_ideal put to 1  !!!'
327                 ELSE
328                    SELECT CASE(isoName(iName))
329                      CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O
330                      CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O
331                      CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O
332                      CASE('HDO');   tnat = tnat_HDO;   alpha_ideal = alpha_ideal_HDO
333                      CASE('HTO');   tnat = tnat_HTO;   alpha_ideal = alpha_ideal_HTO
334                      CASE DEFAULT
335                         CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1)
336                    END SELECT
337                 END IF
338                 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
339              ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
340                 IF(tracers(iq)%iso_iZone == 1) THEN
341                    ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
342                    ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
343                    q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
344                 else !IF(tracers(iq)%iso_iZone == 1) THEN
345                    q(:,:,iq) = 0.
346                 endif !IF(tracers(iq)%iso_iZone == 1) THEN
347              END IF !IF(tracers(iq)%iso_iZone == 0) THEN
348           enddo
349        else
350           q(:,:,:)=0
351        endif ! of if (planet_type=="earth")
352
353        call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
354
355        ! add random perturbation to temperature
356        idum  = -1
357        zz = ran1(idum)
358        idum  = 0
359        do l=1,llm
360           do ij=iip2,ip1jm
361              teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum))
362           enddo
363        enddo
364
365        ! maintain periodicity in longitude
366        do l=1,llm
367           do ij=1,ip1jmp1,iip1
368              teta(ij+iim,l)=teta(ij,l)
369           enddo
370        enddo
371
372     ENDIF ! of IF (.NOT. read_start)
373  endif academic_case
374
375END SUBROUTINE iniacademic
Note: See TracBrowser for help on using the repository browser.