source: LMDZ5/branches/testing/libf/dyn3d_common/inidissip.F90 @ 5425

Last change on this file since 5425 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
RevLine 
[524]1!
[1279]2! $Id: inidissip.F90 2641 2016-09-29 21:26:46Z jyg $
[524]3!
[1489]4SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
[1697]5     tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
[1489]6  !=======================================================================
7  !   initialisation de la dissipation horizontale
8  !=======================================================================
9  !-----------------------------------------------------------------------
10  !   declarations:
11  !   -------------
[524]12
[1502]13  USE control_mod, only : dissip_period,iperiod
[2641]14  USE comconst_mod, ONLY: dissip_deltaz, dissip_factz, dissip_zref, &
15                          dtdiss, dtvr, rad
16  USE comvert_mod, ONLY: preff, presnivs
[1403]17
[1489]18  IMPLICIT NONE
19  include "dimensions.h"
20  include "paramet.h"
21  include "comdissipn.h"
[1490]22  include "iniprint.h"
[524]23
[1490]24  LOGICAL,INTENT(in) :: lstardis
25  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
26  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
27
[1697]28  integer, INTENT(in):: vert_prof_dissip
29  ! Vertical profile of horizontal dissipation
30  ! Allowed values:
[1698]31  ! 0: rational fraction, function of pressure
[1697]32  ! 1: tanh of altitude
33
[1490]34! Local variables:
[1489]35  REAL fact,zvert(llm),zz
[1611]36  REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
37  real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm)
[1489]38  REAL ullm,vllm,umin,vmin,zhmin,zhmax
[1611]39  REAL zllm
[524]40
[1489]41  INTEGER l,ij,idum,ii
42  REAL tetamin
43  REAL pseudoz
[1490]44  character (len=80) :: abort_message
[524]45
[1489]46  REAL ran1
[524]47
48
[1489]49  !-----------------------------------------------------------------------
50  !
51  !   calcul des valeurs propres des operateurs par methode iterrative:
52  !   -----------------------------------------------------------------
[524]53
[1489]54  crot     = -1.
55  cdivu    = -1.
56  cdivh    = -1.
[524]57
[1489]58  !   calcul de la valeur propre de divgrad:
59  !   --------------------------------------
60  idum = 0
61  DO l = 1, llm
62     DO ij = 1, ip1jmp1
[524]63        deltap(ij,l) = 1.
[1489]64     ENDDO
65  ENDDO
[524]66
[1489]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
[524]73
[1489]74  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
[524]75
[1489]76  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
[524]77
[1489]78  IF ( zhmin .GE. zhmax  )     THEN
[1490]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)
[1489]82  ENDIF
[524]83
[1489]84  zllm = ABS( zhmax )
85  DO l = 1,50
86     IF(lstardis) THEN
[1611]87        CALL divgrad2(1,zh,deltap,niterh,divgra)
[1489]88     ELSE
[1611]89        CALL divgrad (1,zh,niterh,divgra)
[1489]90     ENDIF
[524]91
[1611]92     zllm  = ABS(maxval(divgra))
93     zh = divgra / zllm
[1489]94  ENDDO
[524]95
[1489]96  IF(lstardis) THEN
97     cdivh = 1./ zllm
98  ELSE
99     cdivh = zllm ** ( -1./niterh )
100  ENDIF
[524]101
[1489]102  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
103  !   -----------------------------------------------------------------
[1490]104  write(lunout,*)'inidissip: calcul des valeurs propres'
[524]105
[1489]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)
[524]116
[1489]117     CALL minmax(iip1*jjp1,zu,umin,ullm )
118     CALL minmax(iip1*jjm, zv,vmin,vllm )
[524]119
[1489]120     ullm = ABS ( ullm )
121     vllm = ABS ( vllm )
[524]122
[1489]123     DO    l = 1, 50
124        IF(ii.EQ.1) THEN
125           !cccc             CALL covcont( 1,zu,zv,zu,zv )
126           IF(lstardis) THEN
[1611]127              CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy )
[1489]128           ELSE
[1611]129              CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy )
[1489]130           ENDIF
131        ELSE
132           IF(lstardis) THEN
[1611]133              CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy )
[1489]134           ELSE
[1611]135              CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy )
[1489]136           ENDIF
137        ENDIF
[524]138
[1611]139        zllm = max(abs(maxval(gx)), abs(maxval(gy)))
140        zu = gx / zllm
141        zv = gy / zllm
[1489]142     end DO
[524]143
[1489]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
[524]157
[1489]158  end DO
[524]159
[1489]160  !   petit test pour les operateurs non star:
161  !   ----------------------------------------
[524]162
[1489]163  !     IF(.NOT.lstardis) THEN
164  fact    = rad*24./REAL(jjm)
165  fact    = fact*fact
[1490]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
[1489]169  !     ENDIF
[524]170
[1489]171  !-----------------------------------------------------------------------
172  !   variation verticale du coefficient de dissipation:
173  !   --------------------------------------------------
[524]174
[1697]175  if (vert_prof_dissip == 1) then
[1489]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
[524]192
193
[1490]194  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
[524]195
[1489]196  tetamin =  1.e+6
[524]197
[1489]198  DO l=1,llm
199     tetaudiv(l)   = zvert(l)/tetagdiv
200     tetaurot(l)   = zvert(l)/tetagrot
201     tetah(l)      = zvert(l)/tetatemp
[524]202
[1489]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
[524]207
[1502]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
[524]214
[1502]215  dtdiss  = dissip_period * dtvr
216  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
217
[1489]218  DO l = 1,llm
[1490]219     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
[1489]220          dtdiss*tetah(l)
221  ENDDO
[524]222
[1489]223END SUBROUTINE inidissip
Note: See TracBrowser for help on using the repository browser.