[cdftools] 23/228: JMM add cdfpvor-full to the trunk
Alastair McKinstry
mckinstry at moszumanska.debian.org
Fri Jun 12 08:21:24 UTC 2015
This is an automated email from the git hooks/post-receive script.
mckinstry pushed a commit to branch master
in repository cdftools.
commit 4ab218c27d3a39c04b2f8a5b3da33cecc152bb2b
Author: molines <molines at 1055176f-818a-41d9-83e1-73fbe5b947c5>
Date: Mon Apr 19 11:17:58 2010 +0000
JMM add cdfpvor-full to the trunk
git-svn-id: http://servforge.legi.grenoble-inp.fr/svn/CDFTOOLS/trunk@299 1055176f-818a-41d9-83e1-73fbe5b947c5
---
Makefile | 5 +-
cdfpvor-full.f90 | 239 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 243 insertions(+), 1 deletion(-)
diff --git a/Makefile b/Makefile
index 48e8b78..4e500d2 100644
--- a/Makefile
+++ b/Makefile
@@ -18,7 +18,7 @@ EXEC = cdfmoy cdfmoyt cdfmoy_sp cdfstd cdfmoy_sal2_temp2 cdfmoy_annual cdfmoy_
cdfeke cdfrmsssh cdfstdevw cdfstdevts cdflinreg cdfimprovechk\
cdfbn2 cdfbn2-full cdfsig0 cdfsigi cdfsiginsitu cdfbottomsig0 cdfbottomsigi cdfspice\
cdfbottom cdfets cdfcurl cdfw cdfgeo-uv cdfmxl cdfmxl-full\
- cdfrhoproj cdfisopycdep cdfsigintegr cdfpv cdflspv cdfpvor\
+ cdfrhoproj cdfisopycdep cdfsigintegr cdfpv cdflspv cdfpvor cdfpvor-full\
cdfmhst cdfmhst-full cdfvhst cdfvhst-full cdftransportiz cdftransportiz_noheat cdftransportiz-full \
cdftransportizpm \
cdfmasstrp cdfmasstrp-full \
@@ -183,6 +183,9 @@ cdflspv: cdfio.o cdflspv.f90
cdfpvor: cdfio.o cdfpvor.f90
$(F90) cdfpvor.f90 -o cdfpvor cdfio.o eos.o $(FFLAGS)
+cdfpvor-full: cdfio.o cdfpvor-full.f90
+ $(F90) cdfpvor-full.f90 -o cdfpvor-full cdfio.o eos.o $(FFLAGS)
+
cdfkempemekeepe: cdfio.o cdfkempemekeepe.f90
$(F90) cdfkempemekeepe.f90 -o cdfkempemekeepe cdfio.o $(FFLAGS)
diff --git a/cdfpvor-full.f90 b/cdfpvor-full.f90
new file mode 100644
index 0000000..0108b96
--- /dev/null
+++ b/cdfpvor-full.f90
@@ -0,0 +1,239 @@
+PROGRAM cdfpvor_full
+ !!-------------------------------------------------------------------
+ !! *** PROGRAM cdfpvor_full ***
+ !!
+ !! ** Purpose: Compute the Ertel Potential vorticity
+ !! FULL STEP VERSION
+ !!
+ !! ** Method: Try to avoid 3 d arrays : work with 2 levels a a time
+ !! Formula :
+ !! Qpot = drho/dz * ( f + xsi ) = Qstr + Qrel
+ !! * f is the Coriolis factor, computed from the latitudes of the T-grid :
+ !! f(i,j) = 2 * omega * sin ( phit(i,j) * pi / 180 )
+ !!
+ !! * xsi is the relative vorticity (vertical component of the velocity curl),
+ !! computed from the relative vorticity of the F-points interpolated at
+ !! the T-points :
+ !! xsif(i,j) = ( ue(i,j) - ue(i,j+1) - ve(i,j) + ve(i+1,j) ) / areaf(i,j)
+ !! with : ue(i,j) = U(i,j) * e1u(i,j)
+ !! ve(i,j) = V(i,j) * e2v(i,j)
+ !! areaf(i,j) = e1f(i,j) * e2f(i,j)
+ !! xsi(i,j) = ( xsif(i-1,j-1) + xsif(i-1,j) + xsif(i,j-1) + xsif(i,j) ) / 4
+ !! = ( ue(i-1,j-1) + ue(i,j-1) - ue(i-1,j+1) - ue(i,j+1)
+ !! - ve(i-1,j-1) - ve(i-1,j) + ve(i+1,j-1) + ve(i+1,j) )
+ !! / 4 / areat(i,j)
+ !! with : areat(i,j) = e1t(i,j) * e2t(i,j)
+ !!
+ !! units : U, V in m.s-1
+ !! e1u, e2v, e1f, e2f in m
+ !! f, xsi in s-1
+ !! Qpot, Qrel, Qstr in kg.m-4.s-1
+ !!
+ !!
+ !! The brunt-vaisala frequency is computed using the
+ !! polynomial expression of McDougall (1987):
+ !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
+ !! N2 is then insterpolated at T levels
+ !!
+ !! history:
+ !! Original : A.M Treguier december 2005
+ !!--------------------------------------------------------------------
+ !! $Rev: 256 $
+ !! $Date: 2009-07-21 17:49:27 +0200 (Tue, 21 Jul 2009) $
+ !! $Id: cdfpvor.f90 256 2009-07-21 15:49:27Z molines $
+ !!--------------------------------------------------------------
+ !!
+ !! * Modules used
+ USE cdfio
+ USE eos
+
+ !! * Local variables
+ IMPLICIT NONE
+ INTEGER :: jk, jj, ji,jt !: dummy loop index
+ INTEGER :: ierr !: working integer
+ INTEGER :: narg, iargc !:
+ INTEGER :: npiglo,npjglo, npk ,npt !: size of the domain
+ INTEGER :: iup = 1 , idown = 2, itmp
+ INTEGER, DIMENSION(3) :: ipk, id_varout
+ REAL(kind=8) , DIMENSION(:,:), ALLOCATABLE :: e2v, e1u, e1t, e2t, gphit
+ REAL(kind=4) , DIMENSION(:,:), ALLOCATABLE :: un, vn, rotn, zareat, z2fcor , stretch
+ REAL(KIND=4) , DIMENSION (:,:,:), ALLOCATABLE :: ztemp, zsal,zwk !: Array to read 2 layer of data
+ REAL(KIND=4) , DIMENSION (:,:), ALLOCATABLE :: &
+ zn2 , & !: Brunt Vaissala Frequency (N2)
+ tmask, e3w
+ REAL(KIND=4) , DIMENSION (:), ALLOCATABLE :: gdepw ,e3w_0
+ REAL(KIND=4),DIMENSION(:) ,ALLOCATABLE :: tim
+
+ CHARACTER(LEN=256) :: cfilet , cfileu, cfilev, cfileout='pvor.nc' !:
+ CHARACTER(LEN=256) :: coordzgr='mesh_zgr.nc' !:
+ CHARACTER(LEN=256) :: coord ='mesh_hgr.nc' !:
+ TYPE(variable), DIMENSION(3) :: typvar !: structure for attribute
+
+ INTEGER :: ncout
+ INTEGER :: istatus
+ REAL(KIND=4) :: zpi, zomega, rau0sg
+ LOGICAL :: lprint
+
+ rau0sg = 1020/9.81
+ lprint = .false.
+
+ !! Read command line
+ narg= iargc()
+ IF ( narg /= 3 ) THEN
+ PRINT *,' Usage : cdfpvor-full gridT gridU gridV'
+ PRINT *,' FULL STEP VERSION '
+ PRINT *,' Output on pvor.nc, variables vorelvor, vostrvor,vototvor'
+ PRINT *,' Need mesh_zgr.nc and coordinates.nc '
+ STOP
+ ENDIF
+
+ CALL getarg (1, cfilet)
+ CALL getarg (2, cfileu)
+ CALL getarg (3, cfilev)
+
+ npiglo= getdim (cfilet,'x')
+ npjglo= getdim (cfilet,'y')
+ npk = getdim (cfilet,'depth')
+ npt = getdim (cfilet,'time')
+
+ PRINT *, 'npiglo=', npiglo
+ PRINT *, 'npjglo=', npjglo
+ PRINT *, 'npk =', npk
+ PRINT *, 'npt =', npt
+
+ ALLOCATE ( e1u(npiglo,npjglo) , e1t(npiglo,npjglo) )
+ ALLOCATE ( e2v(npiglo,npjglo) , e2t(npiglo,npjglo) )
+ ALLOCATE ( gphit(npiglo,npjglo), z2fcor(npiglo,npjglo))
+ ALLOCATE ( zareat(npiglo,npjglo), stretch(npiglo,npjglo))
+ ALLOCATE ( un(npiglo,npjglo) , vn(npiglo,npjglo) )
+ ALLOCATE ( rotn(npiglo,npjglo) , tmask(npiglo,npjglo) )
+ ALLOCATE (gdepw(npk),tim(npt))
+
+ e1u= getvar(coord, 'e1u', 1,npiglo,npjglo)
+ e1t= getvar(coord, 'e1t', 1,npiglo,npjglo)
+ e2v= getvar(coord, 'e2v', 1,npiglo,npjglo)
+ e2t= getvar(coord, 'e2t', 1,npiglo,npjglo)
+ gphit= getvar(coord, 'gphit', 1,npiglo,npjglo)
+ zpi=ACOS(-1.)
+ zomega = 2*zpi/(3600*24)
+ z2fcor(:,:)=2.0*zomega*SIN(gphit(:,:)*zpi/180.0)
+ zareat(:,:) = 4.*e1t(:,:)*e2t(:,:) ! factor of 4 to normalize relative vorticity
+
+ IF (lprint) print *, ' reading gdepw in file ', trim(coordzgr)
+ gdepw(:) = getvare3(coordzgr, 'gdepw', npk)
+ IF (lprint) print *, ' read gdepw in file ', trim(coordzgr)
+
+ tim=getvar1d(cfileu,'time_counter',1)
+ ierr=putvar1d(ncout,tim,1,'T')
+
+ ALLOCATE (ztemp(npiglo,npjglo,2), zsal(npiglo,npjglo,2))
+ ALLOCATE (zwk(npiglo,npjglo,2) )
+ ALLOCATE (zn2(npiglo,npjglo) , e3w_0(npk), e3w(npiglo,npjglo))
+
+ ! create output fileset
+
+ ipk(:)= npk ! Those three variables are 3D
+ ! define variable name and attribute
+ typvar(1)%name= 'vorelvor'
+ typvar(2)%name= 'vostrvor'
+ typvar(3)%name= 'vototvor'
+ typvar%units='kg.m-4.s-1'
+ typvar%missing_value=0.
+ typvar%valid_min= -1000.
+ typvar%valid_max= 1000.
+ typvar(1)%long_name='Relative_component_of_Ertel_PV'
+ typvar(2)%long_name='Stretching_component_of_Ertel_PV'
+ typvar(3)%long_name='Ertel_potential_vorticity'
+ typvar(1)%short_name='vorelvor'
+ typvar(2)%short_name='vostrvor'
+ typvar(3)%short_name='vototvor'
+ typvar%online_operation='N/A'
+ typvar%axis='TZYX'
+
+
+ ncout =create(cfileout, cfilet, npiglo,npjglo,npk)
+ ierr= createvar (ncout ,typvar,3, ipk,id_varout )
+ ierr= putheadervar(ncout, cfilet,npiglo,npjglo,npk)
+
+
+ DO jt=1,npt
+ ! 2 levels of T and S are required : iup,idown (with respect to W level)
+ ! Compute from bottom to top (for vertical integration)
+ PRINT *,'time=',jt,'(days:',tim(jt)/86400.,')'
+ ztemp(:,:,idown) = getvar(cfilet, 'votemper', npk-1 ,npiglo, npjglo, ktime=jt)
+ zsal( :,:,idown) = getvar(cfilet, 'vosaline', npk-1 ,npiglo, npjglo, ktime=jt)
+ IF (lprint) print *, ' read temperature and salinity at bottom '
+
+ tim=getvar1d(cfilet,'time_counter',npt)
+ ierr=putvar1d(ncout,tim,npt,'T')
+ e3w_0(:) = getvare3(coordzgr, 'e3w',npk)
+
+
+ ! -------------------------------- LOOP OVER LEVELS
+ DO jk = npk-1, 1, -1
+ PRINT *,' level ',jk
+ ! ------------------------------------RELATIVE VORTICITY FIRST
+ IF (lprint) print *, ' trying to read u in file:', trim(cfileu)
+ un(:,:) = getvar(cfileu, 'vozocrtx', jk ,npiglo,npjglo,ktime=jt)
+ IF (lprint) print *, ' trying to read v in file:', trim(cfilev)
+ vn(:,:) = getvar(cfilev, 'vomecrty', jk ,npiglo,npjglo,ktime=jt)
+ un(:,:) = un(:,:)*e1u(:,:) ; vn(:,:) = vn(:,:)*e2v(:,:) ;
+ IF (lprint) print *, ' read u and V OK'
+ ! relative vorticity at T point
+ rotn(:,:) = 0.
+ DO jj = 2, npjglo -1
+ DO ji = 2, npiglo -1
+ rotn(ji,jj) = ( un(ji-1,jj-1) + un(ji,jj-1) &
+ -un(ji-1,jj+1) - un(ji,jj+1) &
+ -vn(ji-1,jj-1) - vn(ji-1,jj) &
+ +vn(ji+1,jj-1) + vn(ji+1,jj)) &
+ /zareat(ji,jj)
+ END DO
+ END DO
+ IF (lprint) print *, ' curl calculated '
+ ! now tmask and Vaisala Frequency bn2
+ IF ( jk > 1) then
+ tmask(:,:)=1.
+ ztemp(:,:,iup)= getvar(cfilet, 'votemper', jk-1 ,npiglo, npjglo,ktime=jt)
+ WHERE(ztemp(:,:,idown) == 0 ) tmask = 0
+ zsal(:,:,iup) = getvar(cfilet, 'vosaline', jk-1 ,npiglo,npjglo,ktime=jt)
+ IF (lprint) print *, ' read temperature and salinity '
+ e3w(:,:) = e3w_0(jk)
+ WHERE (e3w == 0 ) e3w = 1.
+
+ zwk(:,:,iup) = &
+ & eosbn2 ( ztemp,zsal,gdepw(jk),e3w, npiglo,npjglo ,iup,idown)* tmask(:,:)
+ IF (lprint) print *, ' bn2 calculated at w points '
+ !
+ ! now put zn2 at T level (k )
+ WHERE ( zwk(:,:,idown) == 0 )
+ zn2(:,:) = zwk(:,:,iup)
+ ELSEWHERE
+ zn2(:,:) = 0.5 * ( zwk(:,:,iup) + zwk(:,:,idown) ) * tmask(:,:)
+ END WHERE
+ IF (lprint) print *, ' bn2 put back at T points '
+ ENDIF
+ !
+ ! now rotn will be converted to relative vorticity and zn2 to stretching
+ rotn(:,:) = rotn(:,:)* rau0sg * zn2(:,:)
+ stretch(:,:) = zn2(:,:) * rau0sg * z2fcor(:,:)
+
+ ! write the three variables on file at level k
+ ierr = putvar(ncout, id_varout(1) ,rotn*1.e7, jk ,npiglo, npjglo, ktime=jt)
+ ierr = putvar(ncout, id_varout(2) ,stretch*1.e7 , jk, npiglo, npjglo , ktime=jt)
+ ierr = putvar(ncout, id_varout(3) ,(rotn+stretch)*1.e7 , jk, npiglo, npjglo , ktime=jt)
+ IF (lprint) print *, ' three variables written '
+ itmp = idown ; idown = iup ; iup = itmp
+
+ END DO ! loop to next level
+
+ ! set zero at bottom
+ rotn(:,:) = 0
+ ierr = putvar(ncout, id_varout(1) ,rotn, npk ,npiglo, npjglo, ktime=jt)
+ ierr = putvar(ncout, id_varout(2) ,rotn ,npk, npiglo, npjglo , ktime=jt)
+ ierr = putvar(ncout, id_varout(3) ,rotn ,npk, npiglo, npjglo , ktime=jt)
+ END DO ! loop on time
+
+ istatus = closeout(ncout)
+
+END PROGRAM cdfpvor_full
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/cdftools.git
More information about the debian-science-commits
mailing list