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

Timestamp:

06/25/13 09:32:16 (11 years ago)

Author:

dubos

Message:

Towards Lagrangian vertical coordinate (not there yet)

File:

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

codes/icosagcm/trunk/src/disvert.f90

@@ -5 +5 @@
   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
+  ! vertical coordinate : Lagrangian or mass-based
+  INTEGER, SAVE :: caldyn_eta
+  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2
 
 CONTAINS
@@ -91 +95 @@
   END SUBROUTINE init_disvert  
   
+  SUBROUTINE compute_rhodz(comp, ps, rhodz)
+    USE icosa
+    USE omp_para
+    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
+    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
+    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd) :: m, err
+    INTEGER :: l,i,j,ij,dd
+    err=0.
+
+    dd=0
+
+    DO l = ll_begin, ll_end
+       DO j=jj_begin-dd,jj_end+dd
+          DO i=ii_begin-dd,ii_end+dd
+             ij=(j-1)*iim+i
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
+             IF(comp) THEN
+                rhodz(ij,l) = m
+             ELSE
+                err = MAX(err,abs(m-rhodz(ij,l)))
+             END IF
+          ENDDO
+       ENDDO
+    ENDDO
+
+    IF(.NOT. comp) THEN
+       IF(err>1e-10) THEN
+          PRINT *, 'Discrepancy between ps and rhodz detected', err
+          STOP
+       ELSE
+!          PRINT *, 'No discrepancy between ps and rhodz detected'
+       END IF
+    END IF
+
+  END SUBROUTINE compute_rhodz
+
+ 
   SUBROUTINE write_apbp
   USE icosa