source: trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90 @ 270

Last change on this file since 270 was 270, checked in by emillour, 13 years ago

Ehouarn: Mise a jour des dynamiques (seq et ) pour suivre
les changements et developpements de LMDZ5 terrestre
(mise a niveau avec LMDZ5 trunk, rev 1560). Ce qui ne devrais pas changer grand chose au fonctionnement de base du GCM).
Modification importante: correction du bug sur le cumul des flux de masse pour le transport des traceurs (mais dans les faits semble avoir peu d'impact).
(voir "commit_importants.log" pour les details des ajouts et modifications).

File size: 6.9 KB
Line 
1!
2! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $
3!
4SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
5     tetagdiv,tetagrot,tetatemp             )
6  !=======================================================================
7  !   initialisation de la dissipation horizontale
8  !=======================================================================
9  !-----------------------------------------------------------------------
10  !   declarations:
11  !   -------------
12
13  USE control_mod, only : dissip_period,iperiod
14
15  IMPLICIT NONE
16  include "dimensions.h"
17  include "paramet.h"
18  include "comdissipn.h"
19  include "comconst.h"
20  include "comvert.h"
21  include "logic.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! Local variables:
29  REAL fact,zvert(llm),zz
30  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
31  REAL ullm,vllm,umin,vmin,zhmin,zhmax
32  REAL zllm,z1llm
33
34  INTEGER l,ij,idum,ii
35  REAL tetamin
36  REAL Pup
37  character (len=80) :: abort_message
38
39  REAL ran1
40
41
42  !-----------------------------------------------------------------------
43  !
44  !   calcul des valeurs propres des operateurs par methode iterrative:
45  !   -----------------------------------------------------------------
46
47  crot     = -1.
48  cdivu    = -1.
49  cdivh    = -1.
50
51  !   calcul de la valeur propre de divgrad:
52  !   --------------------------------------
53  idum = 0
54  DO l = 1, llm
55     DO ij = 1, ip1jmp1
56        deltap(ij,l) = 1.
57     ENDDO
58  ENDDO
59
60  idum  = -1
61  zh(1) = RAN1(idum)-.5
62  idum  = 0
63  DO ij = 2, ip1jmp1
64     zh(ij) = RAN1(idum) -.5
65  ENDDO
66
67  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
68
69  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
70
71  IF ( zhmin .GE. zhmax  )     THEN
72     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
73     abort_message='probleme generateur alleatoire dans inidissip'
74     call abort_gcm('inidissip',abort_message,1)
75  ENDIF
76
77  zllm = ABS( zhmax )
78  DO l = 1,50
79     IF(lstardis) THEN
80        CALL divgrad2(1,zh,deltap,niterh,zh)
81     ELSE
82        CALL divgrad (1,zh,niterh,zh)
83     ENDIF
84
85     CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
86
87     zllm  = ABS( zhmax )
88     z1llm = 1./zllm
89     DO ij = 1,ip1jmp1
90        zh(ij) = zh(ij)* z1llm
91     ENDDO
92  ENDDO
93
94  IF(lstardis) THEN
95     cdivh = 1./ zllm
96  ELSE
97     cdivh = zllm ** ( -1./niterh )
98  ENDIF
99
100  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
101  !   -----------------------------------------------------------------
102  write(lunout,*)'inidissip: calcul des valeurs propres'
103
104  DO    ii = 1, 2
105     !
106     DO ij = 1, ip1jmp1
107        zu(ij)  = RAN1(idum) -.5
108     ENDDO
109     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
110     DO ij = 1, ip1jm
111        zv(ij) = RAN1(idum) -.5
112     ENDDO
113     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
114
115     CALL minmax(iip1*jjp1,zu,umin,ullm )
116     CALL minmax(iip1*jjm, zv,vmin,vllm )
117
118     ullm = ABS ( ullm )
119     vllm = ABS ( vllm )
120
121     DO    l = 1, 50
122        IF(ii.EQ.1) THEN
123           !cccc             CALL covcont( 1,zu,zv,zu,zv )
124           IF(lstardis) THEN
125              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
126           ELSE
127              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
128           ENDIF
129        ELSE
130           IF(lstardis) THEN
131              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
132           ELSE
133              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
134           ENDIF
135        ENDIF
136
137        CALL minmax(iip1*jjp1,zu,umin,ullm )
138        CALL minmax(iip1*jjm, zv,vmin,vllm )
139
140        ullm = ABS  ( ullm )
141        vllm = ABS  ( vllm )
142
143        zllm  = MAX( ullm,vllm )
144        z1llm = 1./ zllm
145        DO ij = 1, ip1jmp1
146           zu(ij) = zu(ij)* z1llm
147        ENDDO
148        DO ij = 1, ip1jm
149           zv(ij) = zv(ij)* z1llm
150        ENDDO
151     end DO
152
153     IF ( ii.EQ.1 ) THEN
154        IF(lstardis) THEN
155           cdivu  = 1./zllm
156        ELSE
157           cdivu  = zllm **( -1./nitergdiv )
158        ENDIF
159     ELSE
160        IF(lstardis) THEN
161           crot   = 1./ zllm
162        ELSE
163           crot   = zllm **( -1./nitergrot )
164        ENDIF
165     ENDIF
166
167  end DO
168
169  !   petit test pour les operateurs non star:
170  !   ----------------------------------------
171
172  !     IF(.NOT.lstardis) THEN
173  fact    = rad*24./REAL(jjm)
174  fact    = fact*fact
175  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
176  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
177  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
178  !     ENDIF
179
180  !-----------------------------------------------------------------------
181  !   variation verticale du coefficient de dissipation:
182  !   --------------------------------------------------
183
184! First step: going from 1 to dissip_fac_mid (in gcm.def)
185!============
186  DO l=1,llm
187     zz      = 1. - preff/presnivs(l)
188     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
189  ENDDO
190
191  write(lunout,*) 'Dissipation : '
192  write(lunout,*) 'Multiplication de la dissipation en altitude :'
193  write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
194
195! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
196!==========================
197! Utilisation de la fonction tangente hyperbolique pour augmenter
198! arbitrairement la dissipation et donc la stabilite du modele a
199! partir d'une certaine altitude.
200
201!   Le facteur multiplicatif de basse atmosphere etant deja pris
202!   en compte, il faut diviser le facteur multiplicatif de haute
203!   atmosphere par celui-ci.
204
205  if (ok_strato) then
206
207    Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
208    do l=1,llm
209      zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
210                *(1-(0.5*(1+tanh(-6./dissip_deltaz*              &
211               (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
212    enddo
213
214    write(*,*) '  dissip_fac_up =', dissip_fac_up
215    write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
216                   dissip_pupstart,'Pa', dissip_deltaz , '(km)'
217
218  endif
219
220
221  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
222
223  tetamin =  1.e+6
224
225  DO l=1,llm
226     tetaudiv(l)   = zvert(l)/tetagdiv
227     tetaurot(l)   = zvert(l)/tetagrot
228     tetah(l)      = zvert(l)/tetatemp
229
230     IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
231     IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
232     IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
233  ENDDO
234
235  ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
236  IF (dissip_period == 0) THEN
237     dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
238     write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
239     dissip_period = MAX(iperiod,dissip_period)
240  END IF
241
242  dtdiss  = dissip_period * dtvr
243  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
244
245  DO l = 1,llm
246     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
247          dtdiss*tetah(l)
248  ENDDO
249
250END SUBROUTINE inidissip
Note: See TracBrowser for help on using the repository browser.