source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/iniacademic_loc.F90 @ 5736

Last change on this file since 5736 was 5229, checked in by abarral, 14 months ago

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