source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90 @ 5185

Last change on this file since 5185 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
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.1 KB
Line 
1! $Id: inidissip.F90 5159 2024-08-02 19:58:25Z abarral $
2
3SUBROUTINE inidissip(lstardis, nitergdiv, nitergrot, niterh, &
4        tetagdiv, tetagrot, tetatemp, vert_prof_dissip)
5  !=======================================================================
6  !   initialisation de la dissipation horizontale
7  !=======================================================================
8  !-----------------------------------------------------------------------
9  !   declarations:
10  !   -------------
11
12  USE control_mod, ONLY: dissip_period, iperiod
13  USE comconst_mod, ONLY: dissip_deltaz, dissip_factz, dissip_zref, &
14          dtdiss, dtvr, rad
15  USE comvert_mod, ONLY: preff, presnivs
16  USE lmdz_filtreg, ONLY: filtreg
17  USE lmdz_libmath, ONLY: minmax
18  USE lmdz_ran1, ONLY: ran1
19  USE lmdz_iniprint, ONLY: lunout, prt_level
20  USE lmdz_comdissipn, ONLY: tetaudiv, tetaurot, tetah, cdivu, crot, cdivh
21
22USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
23  USE lmdz_paramet
24  IMPLICIT NONE
25
26
27
28  LOGICAL, INTENT(IN) :: lstardis
29  INTEGER, INTENT(IN) :: nitergdiv, nitergrot, niterh
30  REAL, INTENT(IN) :: tetagdiv, tetagrot, tetatemp
31
32  INTEGER, INTENT(IN) :: vert_prof_dissip
33  ! Vertical profile of horizontal dissipation
34  ! Allowed values:
35  ! 0: rational fraction, function of pressure
36  ! 1: tanh of altitude
37
38  ! Local variables:
39  REAL fact, zvert(llm), zz
40  REAL zh(ip1jmp1), zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
41  REAL zv(ip1jm), gy(ip1jm), deltap(ip1jmp1, llm)
42  REAL ullm, vllm, umin, vmin, zhmin, zhmax
43  REAL zllm
44
45  INTEGER l, ij, idum, ii
46  REAL tetamin
47  REAL pseudoz
48  CHARACTER (LEN = 80) :: abort_message
49
50  !-----------------------------------------------------------------------
51
52  !   calcul des valeurs propres des operateurs par methode iterrative:
53  !   -----------------------------------------------------------------
54
55  crot = -1.
56  cdivu = -1.
57  cdivh = -1.
58
59  !   calcul de la valeur propre de divgrad:
60  !   --------------------------------------
61  idum = 0
62  DO l = 1, llm
63    DO ij = 1, ip1jmp1
64      deltap(ij, l) = 1.
65    ENDDO
66  ENDDO
67
68  idum = -1
69  zh(1) = RAN1(idum) - .5
70  idum = 0
71  DO ij = 2, ip1jmp1
72    zh(ij) = RAN1(idum) - .5
73  ENDDO
74
75  CALL filtreg (zh, jjp1, 1, 2, 1, .TRUE., 1)
76
77  CALL minmax(iip1 * jjp1, zh, zhmin, zhmax)
78
79  IF (zhmin >= zhmax)     THEN
80    WRITE(lunout, *)'  Inidissip  zh min max  ', zhmin, zhmax
81    abort_message = 'probleme generateur alleatoire dans inidissip'
82    CALL abort_gcm('inidissip', abort_message, 1)
83  ENDIF
84
85  zllm = ABS(zhmax)
86  DO l = 1, 50
87    IF(lstardis) THEN
88      CALL divgrad2(1, zh, deltap, niterh, divgra)
89    ELSE
90      CALL divgrad (1, zh, niterh, divgra)
91    ENDIF
92
93    zllm = ABS(maxval(divgra))
94    zh = divgra / zllm
95  ENDDO
96
97  IF(lstardis) THEN
98    cdivh = 1. / zllm
99  ELSE
100    cdivh = zllm ** (-1. / niterh)
101  ENDIF
102
103  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
104  !   -----------------------------------------------------------------
105  WRITE(lunout, *)'inidissip: calcul des valeurs propres'
106
107  DO    ii = 1, 2
108
109    DO ij = 1, ip1jmp1
110      zu(ij) = RAN1(idum) - .5
111    ENDDO
112    CALL filtreg (zu, jjp1, 1, 2, 1, .TRUE., 1)
113    DO ij = 1, ip1jm
114      zv(ij) = RAN1(idum) - .5
115    ENDDO
116    CALL filtreg (zv, jjm, 1, 2, 1, .FALSE., 1)
117
118    CALL minmax(iip1 * jjp1, zu, umin, ullm)
119    CALL minmax(iip1 * jjm, zv, vmin, vllm)
120
121    ullm = ABS (ullm)
122    vllm = ABS (vllm)
123
124    DO    l = 1, 50
125      IF(ii==1) THEN
126        !cccc             CALL covcont( 1,zu,zv,zu,zv )
127        IF(lstardis) THEN
128          CALL gradiv2(1, zu, zv, nitergdiv, gx, gy)
129        ELSE
130          CALL gradiv (1, zu, zv, nitergdiv, gx, gy)
131        ENDIF
132      ELSE
133        IF(lstardis) THEN
134          CALL nxgraro2(1, zu, zv, nitergrot, gx, gy)
135        ELSE
136          CALL nxgrarot(1, zu, zv, nitergrot, gx, gy)
137        ENDIF
138      ENDIF
139
140      zllm = max(abs(maxval(gx)), abs(maxval(gy)))
141      zu = gx / zllm
142      zv = gy / zllm
143    end DO
144
145    IF (ii==1) THEN
146      IF(lstardis) THEN
147        cdivu = 1. / zllm
148      ELSE
149        cdivu = zllm **(-1. / nitergdiv)
150      ENDIF
151    ELSE
152      IF(lstardis) THEN
153        crot = 1. / zllm
154      ELSE
155        crot = zllm **(-1. / nitergrot)
156      ENDIF
157    ENDIF
158
159  end DO
160
161  !   petit test pour les operateurs non star:
162  !   ----------------------------------------
163
164  !     IF(.NOT.lstardis) THEN
165  fact = rad * 24. / REAL(jjm)
166  fact = fact * fact
167  WRITE(lunout, *)'inidissip: coef u ', fact / cdivu, 1. / cdivu
168  WRITE(lunout, *)'inidissip: coef r ', fact / crot, 1. / crot
169  WRITE(lunout, *)'inidissip: coef h ', fact / cdivh, 1. / cdivh
170  !     ENDIF
171
172  !-----------------------------------------------------------------------
173  !   variation verticale du coefficient de dissipation:
174  !   --------------------------------------------------
175
176  IF (vert_prof_dissip == 1) THEN
177    DO l = 1, llm
178      pseudoz = 8. * log(preff / presnivs(l))
179      zvert(l) = 1 + &
180              (tanh((pseudoz - dissip_zref) / dissip_deltaz) + 1.) / 2. &
181                      * (dissip_factz - 1.)
182    enddo
183  else
184    DO l = 1, llm
185      zvert(l) = 1.
186    ENDDO
187    fact = 2.
188    DO l = 1, llm
189      zz = 1. - preff / presnivs(l)
190      zvert(l) = fact - (fact - 1.) / (1. + zz * zz)
191    ENDDO
192  ENDIF
193
194  WRITE(lunout, *)'inidissip: Constantes de temps de la diffusion horizontale'
195
196  tetamin = 1.e+6
197
198  DO l = 1, llm
199    tetaudiv(l) = zvert(l) / tetagdiv
200    tetaurot(l) = zvert(l) / tetagrot
201    tetah(l) = zvert(l) / tetatemp
202
203    IF(tetamin> (1. / tetaudiv(l))) tetamin = 1. / tetaudiv(l)
204    IF(tetamin> (1. / tetaurot(l))) tetamin = 1. / tetaurot(l)
205    IF(tetamin> (1. / tetah(l))) tetamin = 1. / tetah(l)
206  ENDDO
207
208  ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
209  IF (dissip_period == 0) THEN
210    dissip_period = INT(tetamin / (2. * dtvr * iperiod)) * iperiod
211    WRITE(lunout, *)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ', tetamin, dtvr, iperiod, dissip_period
212    dissip_period = MAX(iperiod, dissip_period)
213  END IF
214
215  dtdiss = dissip_period * dtvr
216  WRITE(lunout, *)'inidissip: dissip_period=', dissip_period, ' dtdiss=', dtdiss, ' dtvr=', dtvr
217
218  DO l = 1, llm
219    WRITE(lunout, *)zvert(l), dtdiss * tetaudiv(l), dtdiss * tetaurot(l), &
220            dtdiss * tetah(l)
221  ENDDO
222
223END SUBROUTINE inidissip
Note: See TracBrowser for help on using the repository browser.