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

Last change on this file since 5138 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
RevLine 
[524]1!
[1403]2! $Id: iniacademic.F90 5084 2024-07-19 16:40:44Z snguyen $
[524]3!
[1472]4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
[524]5
[2083]6  USE filtreg_mod, ONLY: inifilr
[4325]7  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
[1472]8  USE control_mod, ONLY: day_step,planet_type
[4056]9  use exner_hyb_m, only: exner_hyb
10  use exner_milieu_m, only: exner_milieu
[1454]11#ifdef CPP_IOIPSL
[2083]12  USE IOIPSL, ONLY: getin
[1454]13#else
[1472]14  ! if not using IOIPSL, we still need to use (a local version of) getin
[2083]15  USE ioipsl_getincom, ONLY: getin
[1454]16#endif
[1472]17  USE Write_Field
[2597]18  USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm
[2603]19  USE logic_mod, ONLY: iflag_phys, read_start
[4376]20  USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
[2601]21  USE temps_mod, ONLY: annee_ref, day_ini, day_ref
[2622]22  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
[4120]23  USE readTracFiles_mod, ONLY: addPhase
[5084]24  use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
25  use netcdf, only : NF90_CLOSE, NF90_GET_VAR
[1146]26
[5084]27
[1472]28  !   Author:    Frederic Hourdin      original: 15/01/93
[1474]29  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
30  ! of the American Meteorological Society, 75, 1825.
[524]31
[1472]32  IMPLICIT NONE
[524]33
[1472]34  !   Declararations:
35  !   ---------------
[524]36
[1472]37  include "dimensions.h"
38  include "paramet.h"
39  include "comgeom.h"
40  include "academic.h"
41  include "iniprint.h"
[524]42
[1472]43  !   Arguments:
44  !   ----------
[1146]45
[2083]46  REAL,INTENT(OUT) :: time_0
[1146]47
[2083]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
[524]56
[1472]57  !   Local:
58  !   ------
[524]59
[1472]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
[1474]64  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
65  real tetastrat ! potential temperature in the stratosphere, in K
[1472]66  real tetajl(jjp1,llm)
[4120]67  INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
[524]68
[4419]69  integer :: nid_relief,varid,ierr
70  real, dimension(iip1,jjp1) :: relief
71
[1472]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
[3976]76  REAL phi_pv,dphi_pv,gam_pv,tetanoise   ! Constantes pour polar vortex
[1454]77
[1472]78  real zz,ran1
[4056]79  integer idum
[524]80
[4325]81  REAL zdtvr, tnat, alpha_ideal
[4984]82  LOGICAL,PARAMETER :: tnat1=.true.
[1520]83 
84  character(len=*),parameter :: modname="iniacademic"
85  character(len=80) :: abort_message
[524]86
[2083]87  ! Sanity check: verify that options selected by user are not incompatible
[2087]88  if ((iflag_phys==1).and. .not. read_start) then
[2083]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 
[1472]96  !-----------------------------------------------------------------------
97  ! 1. Initializations for Earth-like case
98  ! --------------------------------------
99  !
100  ! initialize planet radius, rotation rate,...
101  call conf_planete
[524]102
[1472]103  time_0=0.
104  day_ref=1
105  annee_ref=0
[1403]106
[1472]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.
[524]117
[1472]118  if (llm == 1) then
119     ! specific initializations for the shallow water case
120     kappa=1
121  endif
[524]122
[1472]123  CALL iniconst
124  CALL inigeom
125  CALL inifilr
[524]126
[4419]127
128  !------------------------------------------------------------------
[3976]129  ! Initialize pressure and mass field if read_start=.false.
[4419]130  !------------------------------------------------------------------
131
[3976]132  IF (.NOT. read_start) THEN
[4056]133
[4419]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)
[5084]145     if (ierr.EQ.NF90_NOERR) THEN
[4419]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
[3976]168     ! ground geopotential
[4419]169     !phis(:)=0.
170     ps(:)=preff
[3976]171     CALL pression ( ip1jmp1, ap, bp, ps, p       )
[4056]172
[3976]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
[1472]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
[524]189
[1529]190  academic_case: if (iflag_phys >= 2) then
[1472]191     ! initializations
[524]192
[1472]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)
[3976]230     tetanoise=0.005
231     CALL getin('tetanoise',tetanoise)
[524]232
[1472]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
[524]244
[1472]245     ! Potential temperature
246     DO l=1,llm
247        zsig=presnivs(l)/preff
248        tetastrat=ttp*zsig**(-kappa)
249        tetapv=tetastrat
[5084]250        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
[1472]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)
[1505]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
[1472]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
[524]270
[1472]271     !          CALL writefield('theta_eq',tetajl)
[524]272
[1472]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
[3976]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
[1472]291        ! geopotential
292        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
293
[3976]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
[1472]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
[4050]312           do iq=1,nqtot
313              q(:,:,iq)=0.
[4120]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
[2270]316
317              ! CRisi: init des isotopes
318              ! distill de Rayleigh très simplifiée
[4143]319              iName    = tracers(iq)%iso_iName
[4124]320              if (niso <= 0 .OR. iName <= 0) CYCLE
[4056]321              iPhase   = tracers(iq)%iso_iPhase
[4050]322              iqParent = tracers(iq)%iqParent
[4120]323              IF(tracers(iq)%iso_iZone == 0) THEN
[4984]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))) &
[4325]330                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
[4984]331                 endif
[4325]332                 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
[4984]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
[1472]342           enddo
343        else
344           q(:,:,:)=0
345        endif ! of if (planet_type=="earth")
346
[4143]347        call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
[2270]348
[1472]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
[3976]355              teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum))
[1472]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.