source: LMDZ5/branches/IPSLCM6.0.10/libf/dyn3d_common/inidissip.F90 @ 3632

Last change on this file since 3632 was 2641, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2593:2640 into testing branch

  • 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.0 KB
Line 
1!
2! $Id: inidissip.F90 2641 2016-09-29 21:26:46Z acozic $
3!
4SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
5     tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
6  !=======================================================================
7  !   initialisation de la dissipation horizontale
8  !=======================================================================
9  !-----------------------------------------------------------------------
10  !   declarations:
11  !   -------------
12
13  USE control_mod, only : dissip_period,iperiod
14  USE comconst_mod, ONLY: dissip_deltaz, dissip_factz, dissip_zref, &
15                          dtdiss, dtvr, rad
16  USE comvert_mod, ONLY: preff, presnivs
17
18  IMPLICIT NONE
19  include "dimensions.h"
20  include "paramet.h"
21  include "comdissipn.h"
22  include "iniprint.h"
23
24  LOGICAL,INTENT(in) :: lstardis
25  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
26  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
27
28  integer, INTENT(in):: vert_prof_dissip
29  ! Vertical profile of horizontal dissipation
30  ! Allowed values:
31  ! 0: rational fraction, function of pressure
32  ! 1: tanh of altitude
33
34! Local variables:
35  REAL fact,zvert(llm),zz
36  REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
37  real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm)
38  REAL ullm,vllm,umin,vmin,zhmin,zhmax
39  REAL zllm
40
41  INTEGER l,ij,idum,ii
42  REAL tetamin
43  REAL pseudoz
44  character (len=80) :: abort_message
45
46  REAL ran1
47
48
49  !-----------------------------------------------------------------------
50  !
51  !   calcul des valeurs propres des operateurs par methode iterrative:
52  !   -----------------------------------------------------------------
53
54  crot     = -1.
55  cdivu    = -1.
56  cdivh    = -1.
57
58  !   calcul de la valeur propre de divgrad:
59  !   --------------------------------------
60  idum = 0
61  DO l = 1, llm
62     DO ij = 1, ip1jmp1
63        deltap(ij,l) = 1.
64     ENDDO
65  ENDDO
66
67  idum  = -1
68  zh(1) = RAN1(idum)-.5
69  idum  = 0
70  DO ij = 2, ip1jmp1
71     zh(ij) = RAN1(idum) -.5
72  ENDDO
73
74  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
75
76  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
77
78  IF ( zhmin .GE. zhmax  )     THEN
79     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
80     abort_message='probleme generateur alleatoire dans inidissip'
81     call abort_gcm('inidissip',abort_message,1)
82  ENDIF
83
84  zllm = ABS( zhmax )
85  DO l = 1,50
86     IF(lstardis) THEN
87        CALL divgrad2(1,zh,deltap,niterh,divgra)
88     ELSE
89        CALL divgrad (1,zh,niterh,divgra)
90     ENDIF
91
92     zllm  = ABS(maxval(divgra))
93     zh = divgra / zllm
94  ENDDO
95
96  IF(lstardis) THEN
97     cdivh = 1./ zllm
98  ELSE
99     cdivh = zllm ** ( -1./niterh )
100  ENDIF
101
102  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
103  !   -----------------------------------------------------------------
104  write(lunout,*)'inidissip: calcul des valeurs propres'
105
106  DO    ii = 1, 2
107     !
108     DO ij = 1, ip1jmp1
109        zu(ij)  = RAN1(idum) -.5
110     ENDDO
111     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
112     DO ij = 1, ip1jm
113        zv(ij) = RAN1(idum) -.5
114     ENDDO
115     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
116
117     CALL minmax(iip1*jjp1,zu,umin,ullm )
118     CALL minmax(iip1*jjm, zv,vmin,vllm )
119
120     ullm = ABS ( ullm )
121     vllm = ABS ( vllm )
122
123     DO    l = 1, 50
124        IF(ii.EQ.1) THEN
125           !cccc             CALL covcont( 1,zu,zv,zu,zv )
126           IF(lstardis) THEN
127              CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy )
128           ELSE
129              CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy )
130           ENDIF
131        ELSE
132           IF(lstardis) THEN
133              CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy )
134           ELSE
135              CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy )
136           ENDIF
137        ENDIF
138
139        zllm = max(abs(maxval(gx)), abs(maxval(gy)))
140        zu = gx / zllm
141        zv = gy / zllm
142     end DO
143
144     IF ( ii.EQ.1 ) THEN
145        IF(lstardis) THEN
146           cdivu  = 1./zllm
147        ELSE
148           cdivu  = zllm **( -1./nitergdiv )
149        ENDIF
150     ELSE
151        IF(lstardis) THEN
152           crot   = 1./ zllm
153        ELSE
154           crot   = zllm **( -1./nitergrot )
155        ENDIF
156     ENDIF
157
158  end DO
159
160  !   petit test pour les operateurs non star:
161  !   ----------------------------------------
162
163  !     IF(.NOT.lstardis) THEN
164  fact    = rad*24./REAL(jjm)
165  fact    = fact*fact
166  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
167  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
168  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
169  !     ENDIF
170
171  !-----------------------------------------------------------------------
172  !   variation verticale du coefficient de dissipation:
173  !   --------------------------------------------------
174
175  if (vert_prof_dissip == 1) then
176     do l=1,llm
177        pseudoz=8.*log(preff/presnivs(l))
178        zvert(l)=1+ &
179             (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
180             *(dissip_factz-1.)
181     enddo
182  else
183     DO l=1,llm
184        zvert(l)=1.
185     ENDDO
186     fact=2.
187     DO l = 1, llm
188        zz      = 1. - preff/presnivs(l)
189        zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
190     ENDDO
191  endif
192
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.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
204     IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
205     IF( tetamin.GT. (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.