Changeset 156 for codes/icosagcm/trunk/src/disvert.f90

Timestamp:

06/03/13 23:45:04 (11 years ago)

Author:

dubos

Message:

new hydrostatic balance - tested with JW06 - code cleanup to follow

File:

• codes/icosagcm/trunk/src/disvert.f90 (modified) (3 diffs)

codes/icosagcm/trunk/src/disvert.f90

@@ -4 +4 @@
   REAL(rstd), SAVE, POINTER :: bp(:)
   REAL(rstd), SAVE, POINTER :: presnivs(:)
+  REAL(rstd), SAVE, POINTER :: mass_al(:), mass_bl(:), mass_ak(:), mass_bk(:), mass_dak(:), mass_dbk(:)
+  REAL(rstd) :: ptop ! pressure at top of atmosphere l=llm+1
 
 CONTAINS
@@ -17 +19 @@
   IMPLICIT NONE
     CHARACTER(LEN=255) :: disvert_type = 'std'
-    
+    INTEGER :: l
+
     CALL getin("disvert",disvert_type)
     
@@ -61 +64 @@
         
     END SELECT
+
+    ! Convert from pressure coordinate to mass coordinate
+    ! pk = ap(k) + bp(k)*ps(ij) = ptop + g*( mass_ak(k) + mass_bk(k)*ms(ij) )
+    ptop = ap(llm+1)
+    ALLOCATE(mass_al(llm+1))
+    ALLOCATE(mass_bl(llm+1))
+    ALLOCATE(mass_ak(llm))
+    ALLOCATE(mass_bk(llm))
+    ALLOCATE(mass_dak(llm))
+    ALLOCATE(mass_dbk(llm))
+
+    ! FIXME : leave ps for the moment ; change ps to ms later
+    DO l=1,llm+1
+!       mass_al(l) = (ap(l)-ptop)/g
+!       mass_bl(l) = bp(l)/g
+       mass_al(l) = (ap(l)-ptop)
+       mass_bl(l) = bp(l)
+    END DO
+    DO l=1,llm
+       mass_ak(l) = .5*(mass_al(l)+mass_al(l+1))
+       mass_bk(l) = .5*(mass_bl(l)+mass_bl(l+1))
+       mass_dak(l) = mass_al(l)-mass_al(l+1)  ! >0
+       mass_dbk(l) = mass_bl(l)-mass_bl(l+1)  ! >0
+    END DO
 
   END SUBROUTINE init_disvert