[cdftools] 201/228: JMM : bug fix in cdfgeo-uv ; missing allocation for time, and eventual division by 0 at the equator

[Alastair McKinstry]

This is an automated email from the git hooks/post-receive script.

mckinstry pushed a commit to branch master
in repository cdftools.

commit 0036db6f7d9de6354cb4804bdd0546e7f96ce5cb
Author: molines <molines at 1055176f-818a-41d9-83e1-73fbe5b947c5>
Date:   Fri Mar 22 21:14:54 2013 +0000

    JMM : bug fix in cdfgeo-uv ; missing allocation for time, and eventual division by 0 at the equator
    
    
    git-svn-id: http://servforge.legi.grenoble-inp.fr/svn/CDFTOOLS/trunk@649 1055176f-818a-41d9-83e1-73fbe5b947c5
---
 cdfgeo-uv.f90 | 26 ++++++++++++++++++++++----
 1 file changed, 22 insertions(+), 4 deletions(-)

diff --git a/cdfgeo-uv.f90 b/cdfgeo-uv.f90
index 5e192ff..febd6f7 100644
--- a/cdfgeo-uv.f90
+++ b/cdfgeo-uv.f90
@@ -37,6 +37,7 @@ PROGRAM cdfgeo_uv
   INTEGER(KIND=4), DIMENSION(1)             :: id_varoutv     ! varid for vgeo
 
   REAL(KIND=4)                              :: grav           ! gravity
+  REAL(KIND=4)                              :: ffu, ffv       ! coriolis param f at U and V point
   REAL(KIND=4), DIMENSION(:),   ALLOCATABLE :: tim            ! time counter
   REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: e1u, e2v, ff   ! horiz metrics, coriolis (f-point)
   REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: glamu, gphiu   ! longitude latitude u-point
@@ -96,6 +97,11 @@ PROGRAM cdfgeo_uv
   npk    = getdim(cf_tfil, cn_z) 
   npt    = getdim(cf_tfil, cn_t) 
 
+ PRINT *, ' NPIGLO= ', npiglo
+ PRINT *, ' NPJGLO= ', npjglo
+ PRINT *, ' NPK   = ', npk
+ PRINT *, ' NPT   = ', npt
+
   ipk(1)                        = 1
   stypvaru(1)%cname             = TRIM(cn_vozocrtx)
   stypvaru(1)%cunits            = 'm/s'
@@ -119,7 +125,7 @@ PROGRAM cdfgeo_uv
 
   ! Allocate the memory
   ALLOCATE ( e1u(npiglo,npjglo), e2v(npiglo,npjglo) )
-  ALLOCATE ( ff(npiglo,npjglo) )
+  ALLOCATE ( ff(npiglo,npjglo), tim(npt)  )
   ALLOCATE ( glamu(npiglo,npjglo), gphiu(npiglo,npjglo)  )
   ALLOCATE ( glamv(npiglo,npjglo), gphiv(npiglo,npjglo)  )
   ALLOCATE ( un(npiglo,npjglo), vn(npiglo,npjglo)  )
@@ -141,7 +147,7 @@ PROGRAM cdfgeo_uv
   ncoutu = create      (cf_uout, cf_tfil,  npiglo, npjglo, 0                              )
   ierr   = createvar   (ncoutu,  stypvaru, 1,      ipk,    id_varoutu                     )
   ierr   = putheadervar(ncoutu,  cf_tfil,  npiglo, npjglo, 0, pnavlon=glamv, pnavlat=gphiv)
-
+  
   tim  = getvar1d(cf_tfil, cn_vtimec, npt     )
   ierr = putvar1d(ncoutu,  tim,       npt, 'T')
 
@@ -174,10 +180,22 @@ PROGRAM cdfgeo_uv
            END DO
         END DO
         ! e1u and e1v are modified to simplify the computation below
+        ! note that geostrophy is not available near the equator ( f=0)
         DO jj=2, npjglo - 1
            DO ji=2, npiglo - 1
-              e1u(ji,jj)= 2.* grav * umask(ji,jj) / ( ff(ji,jj) + ff(ji,  jj-1) ) / e1u(ji,jj)
-              e2v(ji,jj)= 2.* grav * vmask(ji,jj) / ( ff(ji,jj) + ff(ji-1,jj  ) ) / e2v(ji,jj)
+              ffu = ff(ji,jj) + ff(ji,  jj-1)
+              IF ( ffu /= 0. ) THEN 
+                e1u(ji,jj)= 2.* grav * umask(ji,jj) / ( ffu ) / e1u(ji,jj)
+              ELSE
+                e1u(ji,jj)= 0.  ! spvalue
+              ENDIF
+
+              ffv = ff(ji,jj) + ff(ji-1,jj  )
+              IF ( ffv /= 0. ) THEN 
+                e2v(ji,jj)= 2.* grav * vmask(ji,jj) / ( ffv ) / e2v(ji,jj)
+              ELSE
+                e2v(ji,jj)= 0.  ! spvalue
+              ENDIF
            END DO
         END DO
      END IF

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/cdftools.git