Changeset 1645


Ignore:
Timestamp:
Jul 30, 2012, 6:01:50 PM (12 years ago)
Author:
lguez
Message:

In "disvert", when "vert_sampling" equals "read", read "ap" and "bp"
values instead of "sig". This is more general, it does not assume that :
bp = exp(1 - 1 / sig2)

Tested with gfortran, with "vert_sampling == tropo" and "vert_sampling ==
strato". Identical results.

Location:
LMDZ5/trunk/libf
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d/disvert.F90

    r1623 r1645  
    77  use new_unit_m, only: new_unit
    88  use ioipsl, only: getin
     9  use assert_m, only: assert
    910
    1011  IMPLICIT NONE
    … …  
    2122
    2223  real,intent(in) :: pa, preff
    23   real,intent(out) :: ap(llmp1), bp(llmp1)
     24  real,intent(out) :: ap(llmp1) ! in Pa
     25  real, intent(out):: bp(llmp1)
    2426  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
    2527  real,intent(out) :: presnivs(llm)
    … …  
    7981
    8082     sig(llm+1)=0.
     83
     84     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
     85     bp(llmp1) = 0.
     86
     87     ap = pa * (sig - bp)
    8188  case("tropo")
    8289     DO l = 1, llm
    … …  
    8996        sig(l) = sig(l+1) + dsig(l)
    9097     ENDDO
     98
     99     bp(1)=1.
     100     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     101     bp(llmp1) = 0.
     102
     103     ap(1)=0.
     104     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
    91105  case("strato")
    92106     if (llm==39) then
    … …  
    110124        sig(l) = sig(l+1) + dsig(l)
    111125     ENDDO
     126
     127     bp(1)=1.
     128     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     129     bp(llmp1) = 0.
     130
     131     ap(1)=0.
     132     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
    112133  case("read")
     134     ! Read "ap" and "bp". First line is skipped (title line). "ap"
     135     ! should be in Pa. First couple of values should correspond to
     136     ! the surface, that is : "bp" should be in descending order.
    113137     call new_unit(unit)
    114138     open(unit, file="hybrid.txt", status="old", action="read", &
    … …  
    116140     read(unit, fmt=*) ! skip title line
    117141     do l = 1, llm + 1
    118         read(unit, fmt=*) sig(l)
     142        read(unit, fmt=*) ap(l), bp(l)
    119143     end do
    120144     close(unit)
     145     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
     146          bp(llm + 1) == 0., "disvert: bad ap or bp values")
    121147  case default
    122148     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
    … …  
    130156     nivsig(l)= REAL(l)
    131157  ENDDO
    132 
    133   ! .... Calculs de ap(l) et de bp(l) ....
    134   ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    135 
    136   bp(llmp1) = 0.
    137 
    138   DO l = 1, llm
    139      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    140      ap(l) = pa * ( sig(l) - bp(l) )
    141   ENDDO
    142 
    143   bp(1)=1.
    144   ap(1)=0.
    145 
    146   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    147158
    148159  write(lunout, *)  trim(modname),': BP '

  • LMDZ5/trunk/libf/dyn3dpar/disvert.F90

    r1623 r1645  
    77  use new_unit_m, only: new_unit
    88  use ioipsl, only: getin
     9  use assert_m, only: assert
    910
    1011  IMPLICIT NONE
    … …  
    2122
    2223  real,intent(in) :: pa, preff
    23   real,intent(out) :: ap(llmp1), bp(llmp1)
     24  real,intent(out) :: ap(llmp1) ! in Pa
     25  real, intent(out):: bp(llmp1)
    2426  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
    2527  real,intent(out) :: presnivs(llm)
    … …  
    7981
    8082     sig(llm+1)=0.
     83
     84     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
     85     bp(llmp1) = 0.
     86
     87     ap = pa * (sig - bp)
    8188  case("tropo")
    8289     DO l = 1, llm
    … …  
    8996        sig(l) = sig(l+1) + dsig(l)
    9097     ENDDO
     98
     99     bp(1)=1.
     100     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     101     bp(llmp1) = 0.
     102
     103     ap(1)=0.
     104     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
    91105  case("strato")
    92106     if (llm==39) then
    … …  
    110124        sig(l) = sig(l+1) + dsig(l)
    111125     ENDDO
     126
     127     bp(1)=1.
     128     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     129     bp(llmp1) = 0.
     130
     131     ap(1)=0.
     132     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
    112133  case("read")
     134     ! Read "ap" and "bp". First line is skipped (title line). "ap"
     135     ! should be in Pa. First couple of values should correspond to
     136     ! the surface, that is : "bp" should be in descending order.
    113137     call new_unit(unit)
    114138     open(unit, file="hybrid.txt", status="old", action="read", &
    … …  
    116140     read(unit, fmt=*) ! skip title line
    117141     do l = 1, llm + 1
    118         read(unit, fmt=*) sig(l)
     142        read(unit, fmt=*) ap(l), bp(l)
    119143     end do
    120144     close(unit)
     145     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
     146          bp(llm + 1) == 0., "disvert: bad ap or bp values")
    121147  case default
    122148     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
    … …  
    130156     nivsig(l)= REAL(l)
    131157  ENDDO
    132 
    133   ! .... Calculs de ap(l) et de bp(l) ....
    134   ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    135 
    136   bp(llmp1) = 0.
    137 
    138   DO l = 1, llm
    139      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    140      ap(l) = pa * ( sig(l) - bp(l) )
    141   ENDDO
    142 
    143   bp(1)=1.
    144   ap(1)=0.
    145 
    146   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    147158
    148159  write(lunout, *)  trim(modname),': BP '
Note: See TracChangeset for help on using the changeset viewer.