source: LMDZ6/trunk/libf/dyn3d/iniacademic.f90 @ 5445

Last change on this file since 5445 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
RevLine 
[524]1!
[1403]2! $Id: iniacademic.f90 5292 2024-10-28 15:58:32Z fhourdin $
[524]3!
[1472]4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
[524]5
[5292]6USE academic_mod_h
7    USE iniprint_mod_h
[5281]8  USE comgeom_mod_h
[2083]9  USE filtreg_mod, ONLY: inifilr
[5214]10  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase
[1472]11  USE control_mod, ONLY: day_step,planet_type
[4056]12  use exner_hyb_m, only: exner_hyb
13  use exner_milieu_m, only: exner_milieu
[2083]14  USE IOIPSL, ONLY: getin
[5267]15
[1472]16  USE Write_Field
[2597]17  USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm
[2603]18  USE logic_mod, ONLY: iflag_phys, read_start
[4376]19  USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
[2601]20  USE temps_mod, ONLY: annee_ref, day_ini, day_ref
[2622]21  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
[5084]22  use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
23  use netcdf, only : NF90_CLOSE, NF90_GET_VAR
[5214]24  USE iso_params_mod   ! tnat_* and alpha_ideal_*
[1146]25
[5084]26
[1472]27  !   Author:    Frederic Hourdin      original: 15/01/93
[1474]28  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
29  ! of the American Meteorological Society, 75, 1825.
[524]30
[5271]31  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]32USE paramet_mod_h
[5271]33IMPLICIT NONE
[524]34
[1472]35  !   Declararations:
36  !   ---------------
[524]37
[5271]38
[5272]39
[524]40
[1472]41  !   Arguments:
42  !   ----------
[1146]43
[2083]44  REAL,INTENT(OUT) :: time_0
[1146]45
[2083]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
[524]54
[1472]55  !   Local:
56  !   ------
[524]57
[1472]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
[1474]62  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
63  real tetastrat ! potential temperature in the stratosphere, in K
[1472]64  real tetajl(jjp1,llm)
[4120]65  INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
[524]66
[4419]67  integer :: nid_relief,varid,ierr
68  real, dimension(iip1,jjp1) :: relief
69
[1472]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
[3976]74  REAL phi_pv,dphi_pv,gam_pv,tetanoise   ! Constantes pour polar vortex
[1454]75
[1472]76  real zz,ran1
[4056]77  integer idum
[524]78
[4325]79  REAL zdtvr, tnat, alpha_ideal
[5200]80  LOGICAL :: ltnat1
[1520]81 
82  character(len=*),parameter :: modname="iniacademic"
83  character(len=80) :: abort_message
[524]84
[2083]85  ! Sanity check: verify that options selected by user are not incompatible
[2087]86  if ((iflag_phys==1).and. .not. read_start) then
[2083]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 
[1472]94  !-----------------------------------------------------------------------
95  ! 1. Initializations for Earth-like case
96  ! --------------------------------------
97  !
98  ! initialize planet radius, rotation rate,...
99  call conf_planete
[524]100
[1472]101  time_0=0.
102  day_ref=1
103  annee_ref=0
[1403]104
[1472]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.
[524]115
[1472]116  if (llm == 1) then
117     ! specific initializations for the shallow water case
118     kappa=1
119  endif
[524]120
[1472]121  CALL iniconst
122  CALL inigeom
123  CALL inifilr
[524]124
[4419]125
126  !------------------------------------------------------------------
[3976]127  ! Initialize pressure and mass field if read_start=.false.
[4419]128  !------------------------------------------------------------------
129
[3976]130  IF (.NOT. read_start) THEN
[4056]131
[4419]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)
[5084]143     if (ierr.EQ.NF90_NOERR) THEN
[4419]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
[3976]166     ! ground geopotential
[4419]167     !phis(:)=0.
168     ps(:)=preff
[3976]169     CALL pression ( ip1jmp1, ap, bp, ps, p       )
[4056]170
[3976]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
[1472]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
[524]187
[1529]188  academic_case: if (iflag_phys >= 2) then
[1472]189     ! initializations
[524]190
[1472]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)
[3976]228     tetanoise=0.005
229     CALL getin('tetanoise',tetanoise)
[524]230
[1472]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
[524]242
[1472]243     ! Potential temperature
244     DO l=1,llm
245        zsig=presnivs(l)/preff
246        tetastrat=ttp*zsig**(-kappa)
247        tetapv=tetastrat
[5084]248        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
[1472]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)
[1505]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
[1472]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
[524]268
[1472]269     !          CALL writefield('theta_eq',tetajl)
[524]270
[1472]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
[3976]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
[1472]289        ! geopotential
290        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
291
[3976]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
[1472]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
[5200]309           ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
[1472]310           ! Earth: first two tracers will be water
[4050]311           do iq=1,nqtot
312              q(:,:,iq)=0.
[4120]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
[2270]315
316              ! CRisi: init des isotopes
317              ! distill de Rayleigh très simplifiée
[4143]318              iName    = tracers(iq)%iso_iName
[4124]319              if (niso <= 0 .OR. iName <= 0) CYCLE
[4056]320              iPhase   = tracers(iq)%iso_iPhase
[4050]321              iqParent = tracers(iq)%iqParent
[4120]322              IF(tracers(iq)%iso_iZone == 0) THEN
[5200]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
[5214]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
[5200]337                 END IF
[4325]338                 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
[4984]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
[1472]348           enddo
349        else
350           q(:,:,:)=0
351        endif ! of if (planet_type=="earth")
352
[4143]353        call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
[2270]354
[1472]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
[3976]361              teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum))
[1472]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.