source: LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90 @ 5099

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

Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90
Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad
Create lmdz_yoerad.f90
(lint) Remove "!" on otherwise empty line
  • 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 5099 2024-07-22 19:29:09Z 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,NF90_CLOSE,nf90_get_var
25
26  !   Author:    Frederic Hourdin      original: 15/01/93
27  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
28  ! of the American Meteorological Society, 75, 1825.
29
30  IMPLICIT NONE
31
32  !   Declararations:
33  !   ---------------
34
35  include "dimensions.h"
36  include "paramet.h"
37  include "comgeom.h"
38  include "academic.h"
39  include "iniprint.h"
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,PARAMETER :: tnat1=.true.
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==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<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           ! Earth: first two tracers will be water
310           do iq=1,nqtot
311              q(:,:,iq)=0.
312              IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10
313              IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15
314
315              ! CRisi: init des isotopes
316              ! distill de Rayleigh très simplifiée
317              iName    = tracers(iq)%iso_iName
318              if (niso <= 0 .OR. iName <= 0) CYCLE
319              iPhase   = tracers(iq)%iso_iPhase
320              iqParent = tracers(iq)%iqParent
321              IF(tracers(iq)%iso_iZone == 0) THEN
322                 if (tnat1) then
323                         tnat=1.0
324                         alpha_ideal=1.0
325                         write(*,*) 'Attention dans iniacademic: alpha_ideal=1'
326                 else
327                    IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
328                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
329                 endif
330                 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
331              ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
332                 IF(tracers(iq)%iso_iZone == 1) THEN
333                    ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
334                    ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
335                    q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
336                 else !IF(tracers(iq)%iso_iZone == 1) THEN
337                    q(:,:,iq) = 0.
338                 endif !IF(tracers(iq)%iso_iZone == 1) THEN
339              END IF !IF(tracers(iq)%iso_iZone == 0) THEN
340           enddo
341        else
342           q(:,:,:)=0
343        endif ! of if (planet_type=="earth")
344
345        call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
346
347        ! add random perturbation to temperature
348        idum  = -1
349        zz = ran1(idum)
350        idum  = 0
351        do l=1,llm
352           do ij=iip2,ip1jm
353              teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum))
354           enddo
355        enddo
356
357        ! maintain periodicity in longitude
358        do l=1,llm
359           do ij=1,ip1jmp1,iip1
360              teta(ij+iim,l)=teta(ij,l)
361           enddo
362        enddo
363
364     ENDIF ! of IF (.NOT. read_start)
365  endif academic_case
366
367END SUBROUTINE iniacademic
Note: See TracBrowser for help on using the repository browser.