source: LMDZ6/trunk/libf/dyn3d_common/disvert_noterre.f90 @ 5445

Last change on this file since 5445 was 5285, checked in by abarral, 8 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

  • 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
File size: 9.3 KB
RevLine 
[1520]1! $Id: $
[5246]2SUBROUTINE disvert_noterre
[1520]3
[5246]4  !    Auteur :  F. Forget Y. Wanherdrick, P. Levan
5  !    Nouvelle version 100% Mars !!
6  !    On l'utilise aussi pour Venus et Titan, legerment modifiee.
[1520]7
[5282]8  USE iniprint_mod_h
[5246]9  use IOIPSL
[5267]10
[5246]11  USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt, &
12        nivsig,nivsigs,pa,preff,scaleheight
13  USE comconst_mod, ONLY: kappa
14  USE logic_mod, ONLY: hybrid
[1520]15
[5271]16  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]17USE paramet_mod_h
[5271]18IMPLICIT NONE
[1520]19
[5271]20
[5272]21
[5246]22  !
23  !=======================================================================
[5267]24  !    Discretisation verticale en coordonn�e hybride (ou sigma)
[5246]25  !
26  !=======================================================================
27  !
28  !   declarations:
29  !   -------------
30  !
31  !
32  INTEGER :: l,ll
33  REAL :: snorm
34  REAL :: alpha,beta,gama,delta,deltaz
35  real :: quoi,quand
36  REAL :: zsig(llm),sig(llm+1)
37  INTEGER :: np,ierr
38  integer :: ierr1,ierr2,ierr3,ierr4
39  REAL :: x
[1520]40
[5246]41  REAL :: SSUM
42  EXTERNAL SSUM
43  real :: newsig
44  REAL :: dz0,dz1,nhaut,sig1,esig,csig,zz
45  real :: tt,rr,gg, prevz
46  real :: s(llm),dsig(llm)
[1520]47
[5246]48  integer :: iz
49  real :: z, ps,p
50  character(len=*),parameter :: modname="disvert_noterre"
[1520]51
[5246]52  !
53  !-----------------------------------------------------------------------
54  !
55  ! Initializations:
56  !  pi=2.*ASIN(1.) ! already done in iniconst
[1520]57
[5246]58  hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
59  CALL getin('hybrid',hybrid)
60  write(lunout,*) trim(modname),': hybrid=',hybrid
[1520]61
[5246]62  ! Ouverture possible de fichiers typiquement E.T.
[1520]63
[5246]64     open(99,file="esasig.def",status='old',form='formatted', &
65           iostat=ierr2)
66     if(ierr2.ne.0) then
67          close(99)
68          open(99,file="z2sig.def",status='old',form='formatted', &
69                iostat=ierr4)
70     endif
[1520]71
[5246]72  !-----------------------------------------------------------------------
73  !   cas 1 on lit les options dans esasig.def:
74  !   ----------------------------------------
[1520]75
[5246]76  IF(ierr2.eq.0) then
[1520]77
[5246]78     ! Lecture de esasig.def :
79     ! Systeme peu souple, mais qui respecte en theorie
80     ! La conservation de l'energie (conversion Energie potentielle
81     ! <-> energie cinetique, d'apres la note de Frederic Hourdin...
[1520]82
[5246]83     write(lunout,*)'*****************************'
84     write(lunout,*)'WARNING reading esasig.def'
85     write(lunout,*)'*****************************'
86     READ(99,*) scaleheight
87     READ(99,*) dz0
88     READ(99,*) dz1
89     READ(99,*) nhaut
90     CLOSE(99)
[1520]91
[5246]92     dz0=dz0/scaleheight
93     dz1=dz1/scaleheight
[1520]94
[5246]95     sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
[1520]96
[5246]97     esig=1.
[1520]98
[5246]99     do l=1,20
100        esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.)
101     enddo
102     csig=(1./sig1-1.)/(exp(esig)-1.)
[1520]103
[5246]104     DO L = 2, llm
105        zz=csig*(exp(esig*(l-1.))-1.)
106        sig(l) =1./(1.+zz) &
107              * tanh(.5*(llm+1-l)/nhaut)
108     ENDDO
109     sig(1)=1.
110     sig(llm+1)=0.
111     quoi      = 1. + 2.* kappa
112     s( llm )  = 1.
113     s(llm-1) = quoi
114     IF( llm.gt.2 )  THEN
115        DO  ll = 2, llm-1
116           l         = llm+1 - ll
117           quand     = sig(l+1)/ sig(l)
118           s(l-1)    = quoi * (1.-quand) * s(l)  + quand * s(l+1)
119        ENDDO
120     END IF
121  !
122     snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2)
123     DO l = 1, llm
124        s(l)    = s(l)/ snorm
125     ENDDO
[1520]126
[5246]127  !-----------------------------------------------------------------------
128  !   cas 2 on lit les options dans z2sig.def:
129  !   ----------------------------------------
[1520]130
[5246]131  ELSE IF(ierr4.eq.0) then
132     write(lunout,*)'****************************'
133     write(lunout,*)'Reading z2sig.def'
134     write(lunout,*)'****************************'
[1520]135
[5246]136     READ(99,*) scaleheight
137     do l=1,llm
138        read(99,*) zsig(l)
139     end do
140     CLOSE(99)
[1520]141
[5246]142     sig(1) =1
143     do l=2,llm
144       sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) + &
145             exp(-zsig(l-1)/scaleheight) )
146     end do
147     sig(llm+1) =0
[1520]148
[5246]149  !-----------------------------------------------------------------------
150  ELSE
151     write(lunout,*) 'didn t you forget something ??? '
152     write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
153     stop
154  ENDIF
155  !-----------------------------------------------------------------------
[1520]156
[5246]157  DO l=1,llm
158    nivsigs(l) = REAL(l)
159  ENDDO
[1520]160
[5246]161  DO l=1,llmp1
162    nivsig(l)= REAL(l)
163  ENDDO
[1520]164
165
[5246]166  !-----------------------------------------------------------------------
167  !    ....  Calculs  de ap(l) et de bp(l)  ....
168  !    .........................................
169  !
170  !   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
171  !-----------------------------------------------------------------------
172  !
[1520]173
[5246]174  if (hybrid) then  ! use hybrid coordinates
175     write(lunout,*) "*********************************"
176     write(lunout,*) "Using hybrid vertical coordinates"
177     write(lunout,*)
178     ! Coordonnees hybrides avec mod
179     DO l = 1, llm
[1520]180
[5246]181     call sig_hybrid(sig(l),pa,preff,newsig)
182        bp(l) = EXP( 1. - 1./(newsig**2)  )
183        ap(l) = pa * (newsig - bp(l) )
184     enddo
185     bp(llmp1) = 0.
186     ap(llmp1) = 0.
187  else ! use sigma coordinates
188     write(lunout,*) "********************************"
189     write(lunout,*) "Using sigma vertical coordinates"
190     write(lunout,*)
191     ! Pour ne pas passer en coordonnees hybrides
192     DO l = 1, llm
193        ap(l) = 0.
194        bp(l) = sig(l)
195     ENDDO
196     ap(llmp1) = 0.
197  endif
[1520]198
[5246]199  bp(llmp1) =   0.
[1520]200
[5246]201  write(lunout,*) trim(modname),': BP '
202  write(lunout,*)  bp
203  write(lunout,*) trim(modname),': AP '
204  write(lunout,*)  ap
[1520]205
[5246]206  ! Calcul au milieu des couches :
207  ! WARNING : le choix de placer le milieu des couches au niveau de
[5267]208  ! pression interm�diaire est arbitraire et pourrait etre modifi�.
[5246]209  ! Le calcul du niveau pour la derniere couche
210  ! (on met la meme distance (en log pression)  entre P(llm)
211  ! et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
[5267]212  ! Specifique.  Ce choix est sp�cifi� ici ET dans exner_milieu.F
[1520]213
[5246]214  DO l = 1, llm-1
215   aps(l) =  0.5 *( ap(l) +ap(l+1))
216   bps(l) =  0.5 *( bp(l) +bp(l+1))
217  ENDDO
[1520]218
[5246]219  if (hybrid) then
220     aps(llm) = aps(llm-1)**2 / aps(llm-2)
221     bps(llm) = 0.5*(bp(llm) + bp(llm+1))
222  else
223     bps(llm) = bps(llm-1)**2 / bps(llm-2)
224     aps(llm) = 0.
225  end if
[1520]226
[5246]227  write(lunout,*) trim(modname),': BPs '
228  write(lunout,*)  bps
229  write(lunout,*) trim(modname),': APs'
230  write(lunout,*)  aps
[1520]231
[5246]232  DO l = 1, llm
233   presnivs(l) = aps(l)+bps(l)*preff
234   pseudoalt(l) = -scaleheight*log(presnivs(l)/preff)
235  ENDDO
[1520]236
[5246]237  write(lunout,*)trim(modname),' : PRESNIVS'
238  write(lunout,*)presnivs
239  write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', &
240        'height of ',scaleheight,' km)'
241  write(lunout,*)pseudoalt
[1520]242
[5246]243  ! --------------------------------------------------
244  ! This can be used to plot the vertical discretization
245  ! (> xmgrace -nxy testhybrid.tab )
246  ! --------------------------------------------------
247  ! open (53,file='testhybrid.tab')
248  ! scaleheight=15.5
249  ! do iz=0,34
250  !   z = -5 + min(iz,34-iz)
251  ! approximation of scale height for Venus
252  !   scaleheight = 15.5 - z/55.*10.
253  !   ps = preff*exp(-z/scaleheight)
254  !   zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
255  !   do l=2,llm
256  ! approximation of scale height for Venus
257  !      if (zsig(l-1).le.55.) then
258  !         scaleheight = 15.5 - zsig(l-1)/55.*10.
259  !      else
260  !         scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
261  !      endif
262  !      zsig(l)= zsig(l-1)-scaleheight*
263  !    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
264  !   end do
265  !   write(53,'(I3,50F10.5)') iz, zsig
266  !  end do
267  !  close(53)
268  ! --------------------------------------------------
[1520]269
270
[5246]271  RETURN
272END SUBROUTINE disvert_noterre
[1520]273
[5246]274! ************************************************************
275subroutine sig_hybrid(sig,pa,preff,newsig)
276  ! ----------------------------------------------
277  ! Subroutine utilisee pour calculer des valeurs de sigma modifie
278  ! pour conserver les coordonnees verticales decrites dans
279  ! esasig.def/z2sig.def lors du passage en coordonnees hybrides
280  ! F. Forget 2002
281  ! Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
282  ! L'objectif est de calculer newsig telle que
283  !   (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
[5267]284  ! Cela ne se r�soud pas analytiquement:
285  ! => on r�soud par iterration bourrine
[5246]286  ! ----------------------------------------------
287  ! Information  : where exp(1-1./x**2) become << x
288  !       x      exp(1-1./x**2) /x
289  !       1           1
290  !       0.68       0.5
291  !       0.5        1.E-1
292  !       0.391      1.E-2
293  !       0.333      1.E-3
294  !       0.295      1.E-4
295  !       0.269      1.E-5
296  !       0.248      1.E-6
297  !    => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25
298
299
300  implicit none
301  real :: x1, x2, sig,pa,preff, newsig, F
302  integer :: j
303
304  newsig = sig
305  x1=0
306  x2=1
307  if (sig.ge.1) then
308        newsig= sig
309  else if (sig*preff/pa.ge.0.25) then
310    DO J=1,9999  ! nombre d''iteration max
311      F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
312      ! write(0,*) J, ' newsig =', newsig, ' F= ', F
313      if (F.gt.1) then
314          X2 = newsig
315          newsig=(X1+newsig)*0.5
316      else
317          X1 = newsig
318          newsig=(X2+newsig)*0.5
319      end if
[5267]320      ! Test : on arete lorsque on approxime sig � moins de 0.01 m pr�s
[5246]321      ! (en pseudo altitude) :
322      IF(abs(10.*log(F)).LT.1.E-5) goto 999
323    END DO
324   else   !    if (sig*preff/pa.le.0.25) then
325         newsig= sig*preff/pa
326   end if
[1520]327 999   continue
[5246]328   Return
329END SUBROUTINE sig_hybrid
Note: See TracBrowser for help on using the repository browser.