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

Last change on this file since 5172 was 5159, checked in by abarral, 5 months ago

Put dimensions.h and paramet.h into modules

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