source: LMDZ4/branches/LMDZ4-dev-20091210/libf/dyn3dpar/inidissip.F @ 4999

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

Modif pour la compilation avec fcm sur Vargas SD
Corrections relatives au controle de histins et de l'utilisation des
options INST(X) pour l'appel au routines de IOIPSL FH
Mise en place de diagnostics sur les critères CFL pour l'advection
de traceurs. FH
Controle dans les .def de la dépendance verticale de l'efficacite
de la diffusion. Actif pour le moment uniquement avec ok_strato=y et llm=39. FH
LF

  • 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 1176 2009-06-11 08:54:10Z lebasn $
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.