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

Last change on this file since 5133 was 5128, checked in by abarral, 5 months ago

Correct bug in vlspltqs_loc.f90 from r2270 where we call SSUM with incorrect arguments.
Merge the three different versions of abort_gcm into one
Fix seq, para 3D compilation broken from r5107 onwards
(lint) usual + Remove uneeded fixed-form continuations
  • 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
23
24  !   Author:    Frederic Hourdin      original: 15/01/93
25  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
26  ! of the American Meteorological Society, 75, 1825.
27
28  IMPLICIT NONE
29
30  !   Declararations:
31  !   ---------------
32
33  include "dimensions.h"
34  include "paramet.h"
35  include "comgeom.h"
36  include "academic.h"
37
38  !   Arguments:
39  !   ----------
40
41  REAL, INTENT(OUT) :: time_0
42
43  !   fields
44  REAL, INTENT(OUT) :: vcov(ijb_v:ije_v, llm) ! meridional covariant wind
45  REAL, INTENT(OUT) :: ucov(ijb_u:ije_u, llm) ! zonal covariant wind
46  REAL, INTENT(OUT) :: teta(ijb_u:ije_u, llm) ! potential temperature (K)
47  REAL, INTENT(OUT) :: q(ijb_u:ije_u, llm, nqtot) ! advected tracers (.../kg_of_air)
48  REAL, INTENT(OUT) :: ps(ijb_u:ije_u) ! surface pressure (Pa)
49  REAL, INTENT(OUT) :: masse(ijb_u:ije_u, llm) ! air mass in grid cell (kg)
50  REAL, INTENT(OUT) :: phis(ijb_u:ije_u) ! surface geopotential
51
52  !   Local:
53  !   ------
54
55  REAL, ALLOCATABLE :: vcov_glo(:, :), ucov_glo(:, :), teta_glo(:, :)
56  REAL, ALLOCATABLE :: q_glo(:, :), masse_glo(:, :), ps_glo(:)
57  REAL, ALLOCATABLE :: phis_glo(:)
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
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  ! Initialize pressure and mass field if read_start=.FALSE.
126  IF (.NOT. read_start) THEN
127    ! allocate global fields:
128    !    allocate(vcov_glo(ip1jm,llm))
129
130    allocate(ucov_glo(ip1jmp1, llm))
131    allocate(teta_glo(ip1jmp1, llm))
132    allocate(ps_glo(ip1jmp1))
133    allocate(masse_glo(ip1jmp1, llm))
134    allocate(phis_glo(ip1jmp1))
135
136    ! surface pressure
137    ps_glo(:) = preff
138
139    !------------------------------------------------------------------
140    ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
141    ! du modele.
142    ! On suppose que le fichier relief_in.nc est stoké sur une grille
143    ! iim*jjp1
144    ! Facile a créer à partir de la commande
145    ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
146    !------------------------------------------------------------------
147
148    relief = 0.
149    ierr = nf90_open ('relief_in.nc', nf90_nowrite, nid_relief)
150    IF (ierr==nf90_noerr) THEN
151      ierr = nf90_inq_varid(nid_relief, 'RELIEF', varid)
152      IF (ierr==nf90_noerr) THEN
153        ierr = nf90_get_var(nid_relief, varid, relief(1:iim, 1:jjp1))
154        relief(iip1, :) = relief(1, :)
155      else
156        CALL abort_gcm ('iniacademic', 'variable RELIEF pas la', 1)
157      endif
158    endif
159    ierr = nf90_close (nid_relief)
160
161
162    !------------------------------------------------------------------
163    ! Initialisation du geopotentiel au sol et de la pression
164    !------------------------------------------------------------------
165
166    PRINT*, 'relief=', minval(relief), maxval(relief), 'g=', g
167    do j = 1, jjp1
168      do i = 1, iip1
169        phis_glo((j - 1) * iip1 + i) = g * relief(i, j)
170      enddo
171    enddo
172    PRINT*, 'phis=', minval(phis), maxval(phis), 'g=', g
173
174    CALL pression (ip1jmp1, ap, bp, ps_glo, p)
175    IF (pressure_exner) THEN
176      CALL exner_hyb(ip1jmp1, ps_glo, p, pks, pk)
177    else
178      CALL exner_milieu(ip1jmp1, ps_glo, p, pks, pk)
179    endif
180    CALL massdair(p, masse_glo)
181  ENDIF
182
183  IF (llm == 1) THEN
184    ! initialize fields for the shallow water case, if required
185    IF (.NOT.read_start) THEN
186      phis(ijb_u:ije_u) = 0.
187      q(ijb_u:ije_u, 1:llm, 1:nqtot) = 0
188      CALL sw_case_williamson91_6_loc(vcov, ucov, teta, masse, ps)
189    endif
190  ENDIF
191
192  academic_case: if (iflag_phys >= 2) THEN
193    ! initializations
194
195    ! 1. local parameters
196    ! by convention, winter is in the southern hemisphere
197    ! Geostrophic wind or no wind?
198    ok_geost = .TRUE.
199    CALL getin('ok_geost', ok_geost)
200    ! Constants for Newtonian relaxation and friction
201    k_f = 1.                !friction
202    CALL getin('k_j', k_f)
203    k_f = 1. / (daysec * k_f)
204    k_c_s = 4.  !cooling surface
205    CALL getin('k_c_s', k_c_s)
206    k_c_s = 1. / (daysec * k_c_s)
207    k_c_a = 40. !cooling free atm
208    CALL getin('k_c_a', k_c_a)
209    k_c_a = 1. / (daysec * k_c_a)
210    ! Constants for Teta equilibrium profile
211    teta0 = 315.     ! mean Teta (S.H. 315K)
212    CALL getin('teta0', teta0)
213    ttp = 200.       ! Tropopause temperature (S.H. 200K)
214    CALL getin('ttp', ttp)
215    eps = 0.         ! Deviation to N-S symmetry(~0-20K)
216    CALL getin('eps', eps)
217    delt_y = 60.     ! Merid Temp. Gradient (S.H. 60K)
218    CALL getin('delt_y', delt_y)
219    delt_z = 10.     ! Vertical Gradient (S.H. 10K)
220    CALL getin('delt_z', delt_z)
221    ! Polar vortex
222    ok_pv = .FALSE.
223    CALL getin('ok_pv', ok_pv)
224    phi_pv = -50.            ! Latitude of edge of vortex
225    CALL getin('phi_pv', phi_pv)
226    phi_pv = phi_pv * pi / 180.
227    dphi_pv = 5.             ! Width of the edge
228    CALL getin('dphi_pv', dphi_pv)
229    dphi_pv = dphi_pv * pi / 180.
230    gam_pv = 4.              ! -dT/dz vortex (in K/km)
231    CALL getin('gam_pv', gam_pv)
232    tetanoise = 0.005
233    CALL getin('tetanoise', tetanoise)
234
235    ! 2. Initialize fields towards which to relax
236    ! Friction
237    knewt_g = k_c_a
238    DO l = 1, llm
239      zsig = presnivs(l) / preff
240      knewt_t(l) = (k_c_s - k_c_a) * MAX(0., (zsig - 0.7) / 0.3)
241      kfrict(l) = k_f * MAX(0., (zsig - 0.7) / 0.3)
242    ENDDO
243    DO j = 1, jjp1
244      clat4((j - 1) * iip1 + 1:j * iip1) = cos(rlatu(j))**4
245    ENDDO
246
247    ! Potential temperature
248    DO l = 1, llm
249      zsig = presnivs(l) / preff
250      tetastrat = ttp * zsig**(-kappa)
251      tetapv = tetastrat
252      IF ((ok_pv).AND.(zsig<0.1)) THEN
253        tetapv = tetastrat * (zsig * 10.)**(kappa * cpp * gam_pv / 1000. / g)
254      ENDIF
255      DO j = 1, jjp1
256        ! Troposphere
257        ddsin = sin(rlatu(j))
258        tetajl(j, l) = teta0 - delt_y * ddsin * ddsin + eps * ddsin &
259                - delt_z * (1. - ddsin * ddsin) * log(zsig)
260        IF (planet_type=="giant") THEN
261          tetajl(j, l) = teta0 + (delt_y * &
262                  ((sin(rlatu(j) * 3.14159 * eps + 0.0001))**2)   &
263                  / ((rlatu(j) * 3.14159 * eps + 0.0001)**2))     &
264                  - delt_z * log(zsig)
265        endif
266        ! Profil stratospherique isotherme (+vortex)
267        w_pv = (1. - tanh((rlatu(j) - phi_pv) / dphi_pv)) / 2.
268        tetastrat = tetastrat * (1. - w_pv) + tetapv * w_pv
269        tetajl(j, l) = MAX(tetajl(j, l), tetastrat)
270      ENDDO
271    ENDDO
272
273    !          CALL writefield('theta_eq',tetajl)
274
275    do l = 1, llm
276      do j = 1, jjp1
277        do i = 1, iip1
278          ij = (j - 1) * iip1 + i
279          tetarappel(ij, l) = tetajl(j, l)
280        enddo
281      enddo
282    enddo
283
284    ! 3. Initialize fields (if necessary)
285    IF (.NOT. read_start) THEN
286      ! bulk initialization of temperature
287      IF (iflag_phys>10000) THEN
288        ! Particular case to impose a constant temperature T0=0.01*iflag_phys
289        teta_glo(:, :) = 0.01 * iflag_phys / (pk(:, :) / cpp)
290      ELSE
291        teta_glo(:, :) = tetarappel(:, :)
292      ENDIF
293      ! geopotential
294      CALL geopot(ip1jmp1, teta_glo, pk, pks, phis_glo, phi)
295
296      ! winds
297      IF (ok_geost) THEN
298        CALL ugeostr(phi, ucov_glo)
299      else
300        ucov_glo(:, :) = 0.
301      endif
302      vcov(ijb_v:ije_v, 1:llm) = 0.
303
304      ! bulk initialization of tracers
305      IF (planet_type=="earth") THEN
306        ! Earth: first two tracers will be water
307        do iq = 1, nqtot
308          q(ijb_u:ije_u, :, iq) = 0.
309          IF(tracers(iq)%name == addPhase('H2O', 'g')) q(ijb_u:ije_u, :, iq) = 1.e-10
310          IF(tracers(iq)%name == addPhase('H2O', 'l')) q(ijb_u:ije_u, :, iq) = 1.e-15
311
312          ! CRisi: init des isotopes
313          ! distill de Rayleigh très simplifiée
314          iName = tracers(iq)%iso_iName
315          IF (niso <= 0 .OR. iName <= 0) CYCLE
316          iPhase = tracers(iq)%iso_iPhase
317          iqParent = tracers(iq)%iqParent
318          IF(tracers(iq)%iso_iZone == 0) THEN
319            IF (tnat1) THEN
320              tnat = 1.0
321              alpha_ideal = 1.0
322              WRITE(*, *) 'Attention dans iniacademic: alpha_ideal=1'
323            else
324              IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
325                      CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
326            endif
327            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.)
328          ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
329            IF(tracers(iq)%iso_iZone == 1) THEN ! a verifier.
330              ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
331              ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
332              q(ijb_u:ije_u, :, iq) = q(ijb_u:ije_u, :, iqIsoPha(iName, iPhase))
333            else !IF(tracers(iq)%iso_iZone == 1) THEN
334              q(ijb_u:ije_u, :, iq) = 0.0
335            endif !IF(tracers(iq)%iso_iZone == 1) THEN
336          END IF !IF(tracers(iq)%iso_iZone == 0) THEN
337        enddo
338      else
339        q(ijb_u:ije_u, :, :) = 0
340      endif ! of if (planet_type=="earth")
341
342      CALL check_isotopes(q, ijb_u, ije_u, 'iniacademic_loc')
343
344      ! add random perturbation to temperature
345      idum = -1
346      zz = ran1(idum)
347      idum = 0
348      do l = 1, llm
349        do ij = iip2, ip1jm
350          teta_glo(ij, l) = teta_glo(ij, l) * (1. + tetanoise * ran1(idum))
351        enddo
352      enddo
353
354      ! maintain periodicity in longitude
355      do l = 1, llm
356        do ij = 1, ip1jmp1, iip1
357          teta_glo(ij + iim, l) = teta_glo(ij, l)
358        enddo
359      enddo
360
361      ! copy data from global array to local array:
362      teta(ijb_u:ije_u, :) = teta_glo(ijb_u:ije_u, :)
363      ucov(ijb_u:ije_u, :) = ucov_glo(ijb_u:ije_u, :)
364      !        vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:)
365      masse(ijb_u:ije_u, :) = masse_glo(ijb_u:ije_u, :)
366      ps(ijb_u:ije_u) = ps_glo(ijb_u:ije_u)
367      phis(ijb_u:ije_u) = phis_glo(ijb_u:ije_u)
368
369      deallocate(teta_glo)
370      deallocate(ucov_glo)
371      !        deallocate(vcov_glo)
372      deallocate(masse_glo)
373      deallocate(ps_glo)
374      deallocate(phis_glo)
375    ENDIF ! of IF (.NOT. read_start)
376  ENDIF academic_case
377
378END SUBROUTINE iniacademic_loc
Note: See TracBrowser for help on using the repository browser.