Ignore:
Timestamp:
Dec 23, 2010, 5:38:42 PM (14 years ago)
Author:
lguez
Message:

Conversion to free source form for "disvert" and "iniacademic", no
other change. Bug fix in "fisrtilp": "fraca" could appear in an
expression while not defined.

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90

    r1463 r1472  
    1 !
    21! $Id$
    3 !
    4       SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    52
    6 c    Auteur :  P. Le Van .
    7 c
    8       IMPLICIT NONE
     3SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
    94
    10 #include "dimensions.h"
    11 #include "paramet.h"
    12 #include "iniprint.h"
    13 #include "logic.h"
    14 c
    15 c=======================================================================
    16 c
    17 c
    18 c    s = sigma ** kappa   :  coordonnee  verticale
    19 c    dsig(l)            : epaisseur de la couche l ds la coord.  s
    20 c    sig(l)             : sigma a l'interface des couches l et l-1
    21 c    ds(l)              : distance entre les couches l et l-1 en coord.s
    22 c
    23 c=======================================================================
    24 c
    25       REAL pa,preff
    26       REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
    27       REAL presnivs(llm)
    28 c
    29 c   declarations:
    30 c   -------------
    31 c
    32       REAL sig(llm+1),dsig(llm)
    33        real zzz(1:llm+1)
    34        real dzz(1:llm)
    35       real zk,zkm1,dzk1,dzk2,k0,k1
    36 c
    37       INTEGER l
    38       REAL snorm,dsigmin
    39       REAL alpha,beta,gama,delta,deltaz,h
    40       INTEGER np,ierr
    41       REAL pi,x
     5  ! Auteur : P. Le Van
    426
    43       REAL SSUM
    44 c
    45 c-----------------------------------------------------------------------
    46 c
    47       pi=2.*ASIN(1.)
     7  IMPLICIT NONE
    488
    49       OPEN(99,file='sigma.def',status='old',form='formatted',
    50      s   iostat=ierr)
     9  include "dimensions.h"
     10  include "paramet.h"
     11  include "iniprint.h"
     12  include "logic.h"
    5113
    52 c-----------------------------------------------------------------------
    53 c   cas 1 on lit les options dans sigma.def:
    54 c   ----------------------------------------
     14  ! s = sigma ** kappa : coordonnee verticale
     15  ! dsig(l) : epaisseur de la couche l ds la coord. s
     16  ! sig(l) : sigma a l'interface des couches l et l-1
     17  ! ds(l) : distance entre les couches l et l-1 en coord.s
    5518
    56       IF (ierr.eq.0) THEN
     19  REAL pa, preff
     20  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
     21  REAL presnivs(llm)
    5722
    58       READ(99,*) h           ! hauteur d'echelle 8.
    59       READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
    60       READ(99,*) beta        ! facteur d'acroissement en haut 1.3
    61       READ(99,*) k0          ! nombre de couches dans la transition surf
    62       READ(99,*) k1          ! nombre de couches dans la transition haute
    63       CLOSE(99)
    64       alpha=deltaz/(llm*h)
    65       write(lunout,*)'h,alpha,k0,k1,beta'
     23  REAL sig(llm+1), dsig(llm)
     24  real zk, zkm1, dzk1, dzk2, k0, k1
    6625
    67 c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
     26  INTEGER l
     27  REAL dsigmin
     28  REAL alpha, beta, deltaz, h
     29  INTEGER iostat
     30  REAL pi, x
    6831
    69       alpha=deltaz/tanh(1./k0)*2.
    70       zkm1=0.
    71       sig(1)=1.
    72       do l=1,llm
    73         sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
    74      + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
     32  !-----------------------------------------------------------------------
     33
     34  pi = 2 * ASIN(1.)
     35
     36  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     37
     38  IF (iostat == 0) THEN
     39     ! cas 1 on lit les options dans sigma.def:
     40     READ(99, *) h ! hauteur d'echelle 8.
     41     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
     42     READ(99, *) beta ! facteur d'acroissement en haut 1.3
     43     READ(99, *) k0 ! nombre de couches dans la transition surf
     44     READ(99, *) k1 ! nombre de couches dans la transition haute
     45     CLOSE(99)
     46     alpha=deltaz/(llm*h)
     47     write(lunout, *)'h, alpha, k0, k1, beta'
     48
     49     alpha=deltaz/tanh(1./k0)*2.
     50     zkm1=0.
     51     sig(1)=1.
     52     do l=1, llm
     53        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
     54             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
    7555        zk=-h*log(sig(l+1))
    7656
    7757        dzk1=alpha*tanh(l/k0)
    7858        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
    79         write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
     59        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
    8060        zkm1=zk
    81       enddo
     61     enddo
    8262
    83       sig(llm+1)=0.
     63     sig(llm+1)=0.
    8464
    85 c
    86        DO 2  l = 1, llm
    87        dsig(l) = sig(l)-sig(l+1)
    88    2   CONTINUE
    89 c
     65     DO l = 1, llm
     66        dsig(l) = sig(l)-sig(l+1)
     67     end DO
     68  ELSE
     69     if (ok_strato) then
     70        if (llm==39) then
     71           dsigmin=0.3
     72        else if (llm==50) then
     73           dsigmin=1.
     74        else
     75           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
     76           dsigmin=1.
     77        endif
     78        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
     79     endif
    9080
    91       ELSE
    92 c-----------------------------------------------------------------------
    93 c   cas 2 ancienne discretisation (LMD5...):
    94 c   ----------------------------------------
     81     DO l = 1, llm
     82        x = pi * (l - 0.5) / (llm + 1)
    9583
    96       WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
     84        IF (ok_strato) THEN
     85           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     86                *(1. - tanh(2 * x / pi - 1.))**2 / 4.
     87        ELSE
     88           dsig(l) = 1.0 + 7.0 * SIN(x)**2
     89        ENDIF
     90     ENDDO
     91     dsig = dsig / sum(dsig)
     92     sig(llm+1) = 0.
     93     DO l = llm, 1, -1
     94        sig(l) = sig(l+1) + dsig(l)
     95     ENDDO
     96  ENDIF
    9797
    98       if (ok_strato) then
    99          if (llm==39) then
    100             dsigmin=0.3
    101          else if (llm==50) then
    102             dsigmin=1.
    103          else
    104             WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
    105             dsigmin=1.
    106          endif
    107          WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
    108       endif
     98  DO l=1, llm
     99     nivsigs(l) = REAL(l)
     100  ENDDO
    109101
    110       h=7.
    111       snorm  = 0.
    112       DO l = 1, llm
    113          x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
     102  DO l=1, llmp1
     103     nivsig(l)= REAL(l)
     104  ENDDO
    114105
    115          IF (ok_strato) THEN
    116            dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
    117      &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2       
    118          ELSE
    119            dsig(l) = 1.0 + 7.0 * SIN(x)**2
    120          ENDIF
     106  ! .... Calculs de ap(l) et de bp(l) ....
     107  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    121108
    122          snorm = snorm + dsig(l)
    123       ENDDO
    124       snorm = 1./snorm
    125       DO l = 1, llm
    126          dsig(l) = dsig(l)*snorm
    127       ENDDO
    128       sig(llm+1) = 0.
    129       DO l = llm, 1, -1
    130          sig(l) = sig(l+1) + dsig(l)
    131       ENDDO
     109  bp(llmp1) = 0.
    132110
    133       ENDIF
     111  DO l = 1, llm
     112     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
     113     ap(l) = pa * ( sig(l) - bp(l) )
     114  ENDDO
    134115
     116  bp(1)=1.
     117  ap(1)=0.
    135118
    136       DO l=1,llm
    137         nivsigs(l) = REAL(l)
    138       ENDDO
     119  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    139120
    140       DO l=1,llmp1
    141         nivsig(l)= REAL(l)
    142       ENDDO
     121  write(lunout, *)' BP '
     122  write(lunout, *) bp
     123  write(lunout, *)' AP '
     124  write(lunout, *) ap
    143125
    144 c
    145 c    ....  Calculs  de ap(l) et de bp(l)  ....
    146 c    .........................................
    147 c
    148 c
    149 c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
    150 c
     126  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
     127  write(lunout, *)'couches calcules pour une pression de surface =', preff
     128  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
     129  write(lunout, *)'8km'
     130  DO l = 1, llm
     131     dpres(l) = bp(l) - bp(l+1)
     132     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
     133     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
     134          log(preff/presnivs(l))*8. &
     135          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
     136          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
     137  ENDDO
    151138
    152       bp(llmp1) =   0.
     139  write(lunout, *) 'PRESNIVS '
     140  write(lunout, *) presnivs
    153141
    154       DO l = 1, llm
    155 cc
    156 ccc    ap(l) = 0.
    157 ccc    bp(l) = sig(l)
    158 
    159       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    160       ap(l) = pa * ( sig(l) - bp(l) )
    161 c
    162       ENDDO
    163 
    164       bp(1)=1.
    165       ap(1)=0.
    166 
    167       ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    168 
    169       write(lunout,*)' BP '
    170       write(lunout,*)  bp
    171       write(lunout,*)' AP '
    172       write(lunout,*)  ap
    173 
    174       write(lunout,*)
    175      .'Niveaux de pressions approximatifs aux centres des'
    176       write(lunout,*)'couches calcules pour une pression de surface =',
    177      .                 preff
    178       write(lunout,*)
    179      .     'et altitudes equivalentes pour une hauteur d echelle de'
    180       write(lunout,*)'8km'
    181       DO l = 1, llm
    182        dpres(l) = bp(l) - bp(l+1)
    183        presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
    184        write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
    185      .        log(preff/presnivs(l))*8.
    186      .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
    187      .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
    188       ENDDO
    189 
    190       write(lunout,*)' PRESNIVS '
    191       write(lunout,*)presnivs
    192 
    193       RETURN
    194       END
     142END SUBROUTINE disvert
Note: See TracChangeset for help on using the changeset viewer.