source: LMDZ4/trunk/libf/dyn3d/inidissip.F @ 1351

Last change on this file since 1351 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

  • 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.F 1279 2009-12-10 09:02:56Z musat $
3!
4      SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
5     *                       tetagdiv,tetagrot,tetatemp             )
6c=======================================================================
7c   initialisation de la dissipation horizontale
8c=======================================================================
9c-----------------------------------------------------------------------
10c   declarations:
11c   -------------
12
13      IMPLICIT NONE
14#include "dimensions.h"
15#include "paramet.h"
16#include "comdissipn.h"
17#include "comconst.h"
18#include "comvert.h"
19#include "control.h"
20#include "logic.h"
21
22      LOGICAL lstardis
23      INTEGER nitergdiv,nitergrot,niterh
24      REAL    tetagdiv,tetagrot,tetatemp
25      REAL fact,zvert(llm),zz
26      REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
27      REAL ullm,vllm,umin,vmin,zhmin,zhmax
28      REAL zllm,z1llm
29
30      INTEGER l,ij,idum,ii
31      REAL tetamin
32      REAL pseudoz
33
34      REAL ran1
35
36
37c-----------------------------------------------------------------------
38c
39c   calcul des valeurs propres des operateurs par methode iterrative:
40c   -----------------------------------------------------------------
41
42      crot     = -1.
43      cdivu    = -1.
44      cdivh    = -1.
45
46c   calcul de la valeur propre de divgrad:
47c   --------------------------------------
48      idum = 0
49      DO l = 1, llm
50       DO ij = 1, ip1jmp1
51        deltap(ij,l) = 1.
52       ENDDO
53      ENDDO
54
55      idum  = -1
56      zh(1) = RAN1(idum)-.5
57      idum  = 0
58      DO ij = 2, ip1jmp1
59        zh(ij) = RAN1(idum) -.5
60      ENDDO
61
62      CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
63
64      CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
65
66      IF ( zhmin .GE. zhmax  )     THEN
67         PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
68         STOP'probleme generateur alleatoire dans inidissip'
69      ENDIF
70
71      zllm = ABS( zhmax )
72      DO l = 1,50
73         IF(lstardis) THEN
74            CALL divgrad2(1,zh,deltap,niterh,zh)
75         ELSE
76            CALL divgrad (1,zh,niterh,zh)
77         ENDIF
78
79        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
80
81         zllm  = ABS( zhmax )
82         z1llm = 1./zllm
83         DO ij = 1,ip1jmp1
84            zh(ij) = zh(ij)* z1llm
85         ENDDO
86      ENDDO
87
88      IF(lstardis) THEN
89         cdivh = 1./ zllm
90      ELSE
91         cdivh = zllm ** ( -1./niterh )
92      ENDIF
93
94c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
95c   -----------------------------------------------------------------
96      print*,'calcul des valeurs propres'
97
98      DO  20  ii = 1, 2
99c
100         DO ij = 1, ip1jmp1
101           zu(ij)  = RAN1(idum) -.5
102         ENDDO
103         CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
104         DO ij = 1, ip1jm
105            zv(ij) = RAN1(idum) -.5
106         ENDDO
107         CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
108
109         CALL minmax(iip1*jjp1,zu,umin,ullm )
110         CALL minmax(iip1*jjm, zv,vmin,vllm )
111
112         ullm = ABS ( ullm )
113         vllm = ABS ( vllm )
114
115         DO  5  l = 1, 50
116            IF(ii.EQ.1) THEN
117ccccc             CALL covcont( 1,zu,zv,zu,zv )
118               IF(lstardis) THEN
119                  CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
120               ELSE
121                  CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
122               ENDIF
123            ELSE
124               IF(lstardis) THEN
125                  CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
126               ELSE
127                  CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
128               ENDIF
129            ENDIF
130
131            CALL minmax(iip1*jjp1,zu,umin,ullm )
132            CALL minmax(iip1*jjm, zv,vmin,vllm )
133
134            ullm = ABS  ( ullm )
135            vllm = ABS  ( vllm )
136
137            zllm  = MAX( ullm,vllm )
138            z1llm = 1./ zllm
139            DO ij = 1, ip1jmp1
140              zu(ij) = zu(ij)* z1llm
141            ENDDO
142            DO ij = 1, ip1jm
143               zv(ij) = zv(ij)* z1llm
144            ENDDO
145 5       CONTINUE
146
147         IF ( ii.EQ.1 ) THEN
148            IF(lstardis) THEN
149               cdivu  = 1./zllm
150            ELSE
151               cdivu  = zllm **( -1./nitergdiv )
152            ENDIF
153         ELSE
154            IF(lstardis) THEN
155               crot   = 1./ zllm
156            ELSE
157               crot   = zllm **( -1./nitergrot )
158            ENDIF
159         ENDIF
160
161 20   CONTINUE
162
163c   petit test pour les operateurs non star:
164c   ----------------------------------------
165
166c     IF(.NOT.lstardis) THEN
167         fact    = rad*24./float(jjm)
168         fact    = fact*fact
169         PRINT*,'coef u ', fact/cdivu, 1./cdivu
170         PRINT*,'coef r ', fact/crot , 1./crot
171         PRINT*,'coef h ', fact/cdivh, 1./cdivh
172c     ENDIF
173
174c-----------------------------------------------------------------------
175c   variation verticale du coefficient de dissipation:
176c   --------------------------------------------------
177
178      if (ok_strato .and. llm==39) then
179         do l=1,llm
180            pseudoz=8.*log(preff/presnivs(l))
181            zvert(l)=1+
182     s      (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2.
183     s      *(dissip_factz-1.)
184         enddo
185      else
186         DO l=1,llm
187            zvert(l)=1.
188         ENDDO
189         fact=2.
190         DO l = 1, llm
191            zz      = 1. - preff/presnivs(l)
192            zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
193         ENDDO
194      endif
195
196
197      PRINT*,'Constantes de temps de la diffusion horizontale'
198
199      tetamin =  1.e+6
200
201      DO l=1,llm
202        tetaudiv(l)   = zvert(l)/tetagdiv
203        tetaurot(l)   = zvert(l)/tetagrot
204        tetah(l)      = zvert(l)/tetatemp
205
206        IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
207        IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
208        IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
209      ENDDO
210
211      PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
212      idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
213      PRINT *,' INIDI tetamin idissip ',tetamin,idissip
214      idissip = MAX(iperiod,idissip)
215      dtdiss  = idissip * dtvr
216      PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
217
218      DO l = 1,llm
219         PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
220     *                   dtdiss*tetah(l)
221      ENDDO
222
223c
224      RETURN
225      END
Note: See TracBrowser for help on using the repository browser.