Changeset 1645 for LMDZ5/trunk/libf/dyn3dpar

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 / sig²⁾

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

File:

• LMDZ5/trunk/libf/dyn3dpar/disvert.F90 (modified) (7 diffs)

LMDZ5/trunk/libf/dyn3dpar/disvert.F90

@@ -7 +7 @@
   use new_unit_m, only: new_unit
   use ioipsl, only: getin
+  use assert_m, only: assert
 
   IMPLICIT NONE
@@ -21 +22 @@
 
   real,intent(in) :: pa, preff
-  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: ap(llmp1) ! in Pa
+  real, intent(out):: bp(llmp1)
   real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
   real,intent(out) :: presnivs(llm)
@@ -79 +81 @@
 
      sig(llm+1)=0.
+
+     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
+     bp(llmp1) = 0.
+
+     ap = pa * (sig - bp)
   case("tropo")
      DO l = 1, llm
@@ -89 +96 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
+
+     bp(1)=1.
+     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
+     bp(llmp1) = 0.
+
+     ap(1)=0.
+     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
   case("strato")
      if (llm==39) then
@@ -110 +124 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
+
+     bp(1)=1.
+     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
+     bp(llmp1) = 0.
+
+     ap(1)=0.
+     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
   case("read")
+     ! Read "ap" and "bp". First line is skipped (title line). "ap"
+     ! should be in Pa. First couple of values should correspond to
+     ! the surface, that is : "bp" should be in descending order.
      call new_unit(unit)
      open(unit, file="hybrid.txt", status="old", action="read", &
@@ -116 +140 @@
      read(unit, fmt=*) ! skip title line
      do l = 1, llm + 1
-        read(unit, fmt=*) sig(l)
+        read(unit, fmt=*) ap(l), bp(l)
      end do
      close(unit)
+     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
+          bp(llm + 1) == 0., "disvert: bad ap or bp values")
   case default
      call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
@@ -130 +156 @@
      nivsig(l)= REAL(l)
   ENDDO
-
-  ! .... Calculs de ap(l) et de bp(l) ....
-  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
-
-  bp(llmp1) = 0.
-
-  DO l = 1, llm
-     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
-     ap(l) = pa * ( sig(l) - bp(l) )
-  ENDDO
-
-  bp(1)=1.
-  ap(1)=0.
-
-  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
   write(lunout, *)  trim(modname),': BP '