source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90 @ 5115

Last change on this file since 5115 was 5106, checked in by abarral, 5 months ago

Turn coefils.h into lmdz_coefils.f90
Put filtreg.F90 inside lmdz_filtreg.F90
Turn mod_filtreg_p.F90 into lmdz_filtreg_p.F90
Delete obsolete parafilt.h*
(lint) remove spaces between routine name and args

  • 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 
[5099]1
[1279]2! $Id: inidissip.F90 5106 2024-07-23 20:21:18Z abarral $
[5099]3
[5106]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
[5101]13  USE control_mod, ONLY: dissip_period,iperiod
[2597]14  USE comconst_mod, ONLY: dissip_deltaz, dissip_factz, dissip_zref, &
15                          dtdiss, dtvr, rad
[2600]16  USE comvert_mod, ONLY: preff, presnivs
[5106]17  USE lmdz_filtreg, ONLY: filtreg
[1403]18
[1489]19  IMPLICIT NONE
20  include "dimensions.h"
21  include "paramet.h"
22  include "comdissipn.h"
[1490]23  include "iniprint.h"
[524]24
[1490]25  LOGICAL,INTENT(in) :: lstardis
26  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
27  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
28
[1697]29  integer, INTENT(in):: vert_prof_dissip
30  ! Vertical profile of horizontal dissipation
31  ! Allowed values:
[1698]32  ! 0: rational fraction, function of pressure
[1697]33  ! 1: tanh of altitude
34
[1490]35! Local variables:
[1489]36  REAL fact,zvert(llm),zz
[1611]37  REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
38  real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm)
[1489]39  REAL ullm,vllm,umin,vmin,zhmin,zhmax
[1611]40  REAL zllm
[524]41
[1489]42  INTEGER l,ij,idum,ii
43  REAL tetamin
44  REAL pseudoz
[1490]45  character (len=80) :: abort_message
[524]46
[1489]47  REAL ran1
[524]48
49
[1489]50  !-----------------------------------------------------------------------
[5099]51
[1489]52  !   calcul des valeurs propres des operateurs par methode iterrative:
53  !   -----------------------------------------------------------------
[524]54
[1489]55  crot     = -1.
56  cdivu    = -1.
57  cdivh    = -1.
[524]58
[1489]59  !   calcul de la valeur propre de divgrad:
60  !   --------------------------------------
61  idum = 0
62  DO l = 1, llm
63     DO ij = 1, ip1jmp1
[524]64        deltap(ij,l) = 1.
[1489]65     ENDDO
66  ENDDO
[524]67
[1489]68  idum  = -1
69  zh(1) = RAN1(idum)-.5
70  idum  = 0
71  DO ij = 2, ip1jmp1
72     zh(ij) = RAN1(idum) -.5
73  ENDDO
[524]74
[1489]75  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
[524]76
[1489]77  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
[524]78
[5082]79  IF ( zhmin >= zhmax  )     THEN
[1490]80     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
81     abort_message='probleme generateur alleatoire dans inidissip'
[5101]82     CALL abort_gcm('inidissip',abort_message,1)
[1489]83  ENDIF
[524]84
[1489]85  zllm = ABS( zhmax )
86  DO l = 1,50
87     IF(lstardis) THEN
[1611]88        CALL divgrad2(1,zh,deltap,niterh,divgra)
[1489]89     ELSE
[1611]90        CALL divgrad (1,zh,niterh,divgra)
[1489]91     ENDIF
[524]92
[1611]93     zllm  = ABS(maxval(divgra))
94     zh = divgra / zllm
[1489]95  ENDDO
[524]96
[1489]97  IF(lstardis) THEN
98     cdivh = 1./ zllm
99  ELSE
100     cdivh = zllm ** ( -1./niterh )
101  ENDIF
[524]102
[1489]103  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
104  !   -----------------------------------------------------------------
[1490]105  write(lunout,*)'inidissip: calcul des valeurs propres'
[524]106
[1489]107  DO    ii = 1, 2
[5099]108
[1489]109     DO ij = 1, ip1jmp1
110        zu(ij)  = RAN1(idum) -.5
111     ENDDO
112     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
113     DO ij = 1, ip1jm
114        zv(ij) = RAN1(idum) -.5
115     ENDDO
116     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
[524]117
[1489]118     CALL minmax(iip1*jjp1,zu,umin,ullm )
119     CALL minmax(iip1*jjm, zv,vmin,vllm )
[524]120
[1489]121     ullm = ABS ( ullm )
122     vllm = ABS ( vllm )
[524]123
[1489]124     DO    l = 1, 50
[5082]125        IF(ii==1) THEN
[1489]126           !cccc             CALL covcont( 1,zu,zv,zu,zv )
127           IF(lstardis) THEN
[1611]128              CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy )
[1489]129           ELSE
[1611]130              CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy )
[1489]131           ENDIF
132        ELSE
133           IF(lstardis) THEN
[1611]134              CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy )
[1489]135           ELSE
[1611]136              CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy )
[1489]137           ENDIF
138        ENDIF
[524]139
[1611]140        zllm = max(abs(maxval(gx)), abs(maxval(gy)))
141        zu = gx / zllm
142        zv = gy / zllm
[1489]143     end DO
[524]144
[5082]145     IF ( ii==1 ) THEN
[1489]146        IF(lstardis) THEN
147           cdivu  = 1./zllm
148        ELSE
149           cdivu  = zllm **( -1./nitergdiv )
150        ENDIF
151     ELSE
152        IF(lstardis) THEN
153           crot   = 1./ zllm
154        ELSE
155           crot   = zllm **( -1./nitergrot )
156        ENDIF
157     ENDIF
[524]158
[1489]159  end DO
[524]160
[1489]161  !   petit test pour les operateurs non star:
162  !   ----------------------------------------
[524]163
[1489]164  !     IF(.NOT.lstardis) THEN
165  fact    = rad*24./REAL(jjm)
166  fact    = fact*fact
[1490]167  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
168  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
169  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
[1489]170  !     ENDIF
[524]171
[1489]172  !-----------------------------------------------------------------------
173  !   variation verticale du coefficient de dissipation:
174  !   --------------------------------------------------
[524]175
[1697]176  if (vert_prof_dissip == 1) then
[1489]177     do l=1,llm
178        pseudoz=8.*log(preff/presnivs(l))
179        zvert(l)=1+ &
180             (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
181             *(dissip_factz-1.)
182     enddo
183  else
184     DO l=1,llm
185        zvert(l)=1.
186     ENDDO
187     fact=2.
188     DO l = 1, llm
189        zz      = 1. - preff/presnivs(l)
190        zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
191     ENDDO
192  endif
[524]193
194
[1490]195  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
[524]196
[1489]197  tetamin =  1.e+6
[524]198
[1489]199  DO l=1,llm
200     tetaudiv(l)   = zvert(l)/tetagdiv
201     tetaurot(l)   = zvert(l)/tetagrot
202     tetah(l)      = zvert(l)/tetatemp
[524]203
[5082]204     IF( tetamin> (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
205     IF( tetamin> (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
206     IF( tetamin> (1./   tetah(l)) ) tetamin = 1./    tetah(l)
[1489]207  ENDDO
[524]208
[1502]209  ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
210  IF (dissip_period == 0) THEN
211     dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
212     write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
213     dissip_period = MAX(iperiod,dissip_period)
214  END IF
[524]215
[1502]216  dtdiss  = dissip_period * dtvr
217  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
218
[1489]219  DO l = 1,llm
[1490]220     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
[1489]221          dtdiss*tetah(l)
222  ENDDO
[524]223
[1489]224END SUBROUTINE inidissip
Note: See TracBrowser for help on using the repository browser.