Ignore:
Timestamp:
Oct 21, 2024, 2:58:45 PM (23 hours ago)
Author:
abarral
Message:

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d_common/disvert_noterre.F90

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