[Dart-dev] [5743] DART/branches/development: This is a uniform PV two-surface QG+1 spectral model contributed
nancy at ucar.edu
nancy at ucar.edu
Thu May 31 16:47:55 MDT 2012
Revision: 5743
Author: nancy
Date: 2012-05-31 16:47:55 -0600 (Thu, 31 May 2012)
Log Message:
-----------
This is a uniform PV two-surface QG+1 spectral model contributed
by Rahul Mahajan, based on a model from Greg Hakim. see the
model_description file for an important note about a change
needed in the obs_model_mod.f90 code. rahul also contributed
updated matlab diagnostic files with cases for the 'sqg' model,
but i haven't merged them yet with the most recent changes
tim made in the diags. i'll commit those as soon as i can.
Added Paths:
-----------
DART/branches/development/models/sqg/
DART/branches/development/models/sqg/fft.f90
DART/branches/development/models/sqg/model_description
DART/branches/development/models/sqg/model_mod.f90
DART/branches/development/models/sqg/spectral_mod.f90
DART/branches/development/models/sqg/sqg.f90
DART/branches/development/models/sqg/sqg_mod.f90
DART/branches/development/models/sqg/work/
DART/branches/development/models/sqg/work/input.nml
DART/branches/development/models/sqg/work/mkmf_closest_member_tool
DART/branches/development/models/sqg/work/mkmf_create_fixed_network_seq
DART/branches/development/models/sqg/work/mkmf_create_obs_sequence
DART/branches/development/models/sqg/work/mkmf_filter
DART/branches/development/models/sqg/work/mkmf_integrate_model
DART/branches/development/models/sqg/work/mkmf_obs_diag
DART/branches/development/models/sqg/work/mkmf_obs_seq_to_netcdf
DART/branches/development/models/sqg/work/mkmf_obs_sequence_tool
DART/branches/development/models/sqg/work/mkmf_perfect_model_obs
DART/branches/development/models/sqg/work/mkmf_preprocess
DART/branches/development/models/sqg/work/mkmf_restart_file_tool
DART/branches/development/models/sqg/work/mkmf_sqg
DART/branches/development/models/sqg/work/mkmf_wakeup_filter
DART/branches/development/models/sqg/work/mkmf_write_sqg_restart
DART/branches/development/models/sqg/work/obs_network.py
DART/branches/development/models/sqg/work/obs_seq.in
DART/branches/development/models/sqg/work/obs_seq.out
DART/branches/development/models/sqg/work/path_names_closest_member_tool
DART/branches/development/models/sqg/work/path_names_create_fixed_network_seq
DART/branches/development/models/sqg/work/path_names_create_obs_sequence
DART/branches/development/models/sqg/work/path_names_filter
DART/branches/development/models/sqg/work/path_names_integrate_model
DART/branches/development/models/sqg/work/path_names_obs_diag
DART/branches/development/models/sqg/work/path_names_obs_seq_to_netcdf
DART/branches/development/models/sqg/work/path_names_obs_sequence_tool
DART/branches/development/models/sqg/work/path_names_perfect_model_obs
DART/branches/development/models/sqg/work/path_names_preprocess
DART/branches/development/models/sqg/work/path_names_restart_file_tool
DART/branches/development/models/sqg/work/path_names_sqg
DART/branches/development/models/sqg/work/path_names_wakeup_filter
DART/branches/development/models/sqg/work/path_names_write_sqg_restart
DART/branches/development/models/sqg/work/quickbuild.csh
DART/branches/development/models/sqg/work/sqgRestart.nc
DART/branches/development/models/sqg/write_sqg_restart.f90
DART/branches/development/obs_def/obs_def_sqg_mod.f90
-------------- next part --------------
Added: DART/branches/development/models/sqg/fft.f90
===================================================================
--- DART/branches/development/models/sqg/fft.f90 (rev 0)
+++ DART/branches/development/models/sqg/fft.f90 2012-05-31 22:47:55 UTC (rev 5743)
@@ -0,0 +1,563 @@
+! ===========================================================
+! <next few lines under version control, D O N O T E D I T>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+! ===========================================================
+
+! ======================================================================
+! NIST Guide to Available Math Software.
+! Source for module FFT from package GO.
+! Retrieved from NETLIB on Wed Jul 5 11:50:07 1995.
+! ======================================================================
+
+ subroutine fft(a,b,ntot,n,nspan,isn)
+
+! multivariate complex fourier transform, computed in place
+! using mixed-radix fast fourier transform algorithm.
+! by r. c. singleton, stanford research institute, sept. 1968
+! arrays a and b originally hold the real and imaginary
+! components of the data, and return the real and
+! imaginary components of the resulting fourier coefficients.
+! multivariate data is indexed according to the fortran
+! array element successor function, without limit
+! on the number of implied multiple subscripts.
+! the subroutine is called once for each variate.
+! the calls for a multivariate transform may be in any order.
+! ntot is the total number of complex data values.
+! n is the dimension of the current variable.
+! nspan/n is the spacing of consecutive data values
+! while indexing the current variable.
+! the sign of isn determines the sign of the complex
+! exponential, and the magnitude of isn is normally one.
+! a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
+! is computed by
+! call fft(a,b,n1*n2*n3,n1,n1,1)
+! call fft(a,b,n1*n2*n3,n2,n1*n2,1)
+! call fft(a,b,n1*n2*n3,n3,n1*n2*n3,1)
+! for a single-variate transform,
+! ntot = n = nspan = (number of complex data values), e.g.
+! call fft(a,b,n,n,n,1)
+! the data can alternatively be stored in a single complex array c
+! in standard fortran fashion, i.e. alternating real and imaginary
+! parts. then with most fortran compilers, the complex array c can
+! be equivalenced to a real array a, the magnitude of isn changed
+! to two to give correct indexing increment, and a(1) and a(2) used
+! to pass the initial addresses for the sequences of real and
+! imaginary values, e.g.
+! complex c(ntot)
+! real a(2*ntot)
+! equivalence (c(1),a(1))
+! call fft(a(1),a(2),ntot,n,nspan,2)
+! arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
+! are used for temporary storage. if the available storage
+! is insufficient, the program is terminated by a stop.
+! maxf must be .ge. the maximum prime factor of n.
+! maxp must be .gt. the number of prime factors of n.
+! in addition, if the square-free portion k of n has two or
+! more prime factors, then maxp must be .ge. k-1.
+! dimension a(1),b(1)
+ dimension a(*),b(*)
+! array storage in nfac for a maximum of 15 prime factors of n.
+! if n has more than one square-free factor, the product of the
+! square-free factors must be .le. 210
+ dimension nfac(101),np(209)
+! array storage for maximum prime factor of 23
+ dimension at(53),ck(53),bt(53),sk(53)
+ equivalence (i,ii)
+! the following two constants should agree with the array dimensions.
+ maxp=209
+ if(n .lt. 2) return
+ inc=isn
+ c72=0.30901699437494742
+ s72=0.95105651629515357
+ s120=0.86602540378443865
+ rad=6.2831853071796
+ if(isn .ge. 0) go to 10
+ s72=-s72
+ s120=-s120
+ rad=-rad
+ inc=-inc
+ 10 nt=inc*ntot
+ ks=inc*nspan
+ kspan=ks
+ nn=nt-inc
+ jc=ks/n
+ radf=rad*float(jc)*0.5
+ i=0
+ jf=0
+! determine the factors of n
+ m=0
+ k=n
+ go to 20
+ 15 m=m+1
+ nfac(m)=4
+ k=k/16
+ 20 if(k-(k/16)*16 .eq. 0) go to 15
+ j=3
+ jj=9
+ go to 30
+ 25 m=m+1
+ nfac(m)=j
+ k=k/jj
+ 30 if(mod(k,jj) .eq. 0) go to 25
+ j=j+2
+ jj=j**2
+ if(jj .le. k) go to 30
+ if(k .gt. 4) go to 40
+ kt=m
+ nfac(m+1)=k
+ if(k .ne. 1) m=m+1
+ go to 80
+ 40 if(k-(k/4)*4 .ne. 0) go to 50
+ m=m+1
+ nfac(m)=2
+ k=k/4
+ 50 kt=m
+ j=2
+ 60 if(mod(k,j) .ne. 0) go to 70
+ m=m+1
+ nfac(m)=j
+ k=k/j
+ 70 j=((j+1)/2)*2+1
+ if(j .le. k) go to 60
+ 80 if(kt .eq. 0) go to 100
+ j=kt
+ 90 m=m+1
+ nfac(m)=nfac(j)
+ j=j-1
+ if(j .ne. 0) go to 90
+! compute fourier transform
+ 100 sd=radf/float(kspan)
+ cd=2.0*sin(sd)**2
+ sd=sin(sd+sd)
+ kk=1
+ i=i+1
+ if(nfac(i) .ne. 2) go to 400
+! transform for factor of 2 (including rotation factor)
+ kspan=kspan/2
+ k1=kspan+2
+ 210 k2=kk+kspan
+ ak=a(k2)
+ bk=b(k2)
+ a(k2)=a(kk)-ak
+ b(k2)=b(kk)-bk
+ a(kk)=a(kk)+ak
+ b(kk)=b(kk)+bk
+ kk=k2+kspan
+ if(kk .le. nn) go to 210
+ kk=kk-nn
+ if(kk .le. jc) go to 210
+ if(kk .gt. kspan) go to 800
+ 220 c1=1.0-cd
+ s1=sd
+ 230 k2=kk+kspan
+ ak=a(kk)-a(k2)
+ bk=b(kk)-b(k2)
+ a(kk)=a(kk)+a(k2)
+ b(kk)=b(kk)+b(k2)
+ a(k2)=c1*ak-s1*bk
+ b(k2)=s1*ak+c1*bk
+ kk=k2+kspan
+ if(kk .lt. nt) go to 230
+ k2=kk-nt
+ c1=-c1
+ kk=k1-k2
+ if(kk .gt. k2) go to 230
+ ak=c1-(cd*c1+sd*s1)
+ s1=(sd*c1-cd*s1)+s1
+ c1=2.0-(ak**2+s1**2)
+ s1=c1*s1
+ c1=c1*ak
+ kk=kk+jc
+ if(kk .lt. k2) go to 230
+ k1=k1+inc+inc
+ kk=(k1-kspan)/2+jc
+ if(kk .le. jc+jc) go to 220
+ go to 100
+! transform for factor of 3 (optional code)
+ 320 k1=kk+kspan
+ k2=k1+kspan
+ ak=a(kk)
+ bk=b(kk)
+ aj=a(k1)+a(k2)
+ bj=b(k1)+b(k2)
+ a(kk)=ak+aj
+ b(kk)=bk+bj
+ ak=-0.5*aj+ak
+ bk=-0.5*bj+bk
+ aj=(a(k1)-a(k2))*s120
+ bj=(b(k1)-b(k2))*s120
+ a(k1)=ak-bj
+ b(k1)=bk+aj
+ a(k2)=ak+bj
+ b(k2)=bk-aj
+ kk=k2+kspan
+ if(kk .lt. nn) go to 320
+ kk=kk-nn
+ if(kk .le. kspan) go to 320
+ go to 700
+! transform for factor of 4
+ 400 if(nfac(i) .ne. 4) go to 600
+ kspnn=kspan
+ kspan=kspan/4
+ 410 c1=1.0
+ s1=0
+ 420 k1=kk+kspan
+ k2=k1+kspan
+ k3=k2+kspan
+ akp=a(kk)+a(k2)
+ akm=a(kk)-a(k2)
+ ajp=a(k1)+a(k3)
+ ajm=a(k1)-a(k3)
+ a(kk)=akp+ajp
+ ajp=akp-ajp
+ bkp=b(kk)+b(k2)
+ bkm=b(kk)-b(k2)
+ bjp=b(k1)+b(k3)
+ bjm=b(k1)-b(k3)
+ b(kk)=bkp+bjp
+ bjp=bkp-bjp
+ if(isn .lt. 0) go to 450
+ akp=akm-bjm
+ akm=akm+bjm
+ bkp=bkm+ajm
+ bkm=bkm-ajm
+ if(s1 .eq. 0) go to 460
+ 430 a(k1)=akp*c1-bkp*s1
+ b(k1)=akp*s1+bkp*c1
+ a(k2)=ajp*c2-bjp*s2
+ b(k2)=ajp*s2+bjp*c2
+ a(k3)=akm*c3-bkm*s3
+ b(k3)=akm*s3+bkm*c3
+ kk=k3+kspan
+ if(kk .le. nt) go to 420
+ 440 c2=c1-(cd*c1+sd*s1)
+ s1=(sd*c1-cd*s1)+s1
+ c1=2.0-(c2**2+s1**2)
+ s1=c1*s1
+ c1=c1*c2
+ c2=c1**2-s1**2
+ s2=2.0*c1*s1
+ c3=c2*c1-s2*s1
+ s3=c2*s1+s2*c1
+ kk=kk-nt+jc
+ if(kk .le. kspan) go to 420
+ kk=kk-kspan+inc
+ if(kk .le. jc) go to 410
+ if(kspan .eq. jc) go to 800
+ go to 100
+ 450 akp=akm+bjm
+ akm=akm-bjm
+ bkp=bkm-ajm
+ bkm=bkm+ajm
+ if(s1 .ne. 0) go to 430
+ 460 a(k1)=akp
+ b(k1)=bkp
+ a(k2)=ajp
+ b(k2)=bjp
+ a(k3)=akm
+ b(k3)=bkm
+ kk=k3+kspan
+ if(kk .le. nt) go to 420
+ go to 440
+! transform for factor of 5 (optional code)
+ 510 c2=c72**2-s72**2
+ s2=2.0*c72*s72
+ 520 k1=kk+kspan
+ k2=k1+kspan
+ k3=k2+kspan
+ k4=k3+kspan
+ akp=a(k1)+a(k4)
+ akm=a(k1)-a(k4)
+ bkp=b(k1)+b(k4)
+ bkm=b(k1)-b(k4)
+ ajp=a(k2)+a(k3)
+ ajm=a(k2)-a(k3)
+ bjp=b(k2)+b(k3)
+ bjm=b(k2)-b(k3)
+ aa=a(kk)
+ bb=b(kk)
+ a(kk)=aa+akp+ajp
+ b(kk)=bb+bkp+bjp
+ ak=akp*c72+ajp*c2+aa
+ bk=bkp*c72+bjp*c2+bb
+ aj=akm*s72+ajm*s2
+ bj=bkm*s72+bjm*s2
+ a(k1)=ak-bj
+ a(k4)=ak+bj
+ b(k1)=bk+aj
+ b(k4)=bk-aj
+ ak=akp*c2+ajp*c72+aa
+ bk=bkp*c2+bjp*c72+bb
+ aj=akm*s2-ajm*s72
+ bj=bkm*s2-bjm*s72
+ a(k2)=ak-bj
+ a(k3)=ak+bj
+ b(k2)=bk+aj
+ b(k3)=bk-aj
+ kk=k4+kspan
+ if(kk .lt. nn) go to 520
+ kk=kk-nn
+ if(kk .le. kspan) go to 520
+ go to 700
+! transform for odd factors
+ 600 k=nfac(i)
+ kspnn=kspan
+ kspan=kspan/k
+ if(k .eq. 3) go to 320
+ if(k .eq. 5) go to 510
+ if(k .eq. jf) go to 640
+ jf=k
+ s1=rad/float(k)
+ c1=cos(s1)
+ s1=sin(s1)
+ if(jf .gt. maxf) go to 998
+ ck(jf)=1.0
+ sk(jf)=0.0
+ j=1
+ 630 ck(j)=ck(k)*c1+sk(k)*s1
+ sk(j)=ck(k)*s1-sk(k)*c1
+ k=k-1
+ ck(k)=ck(j)
+ sk(k)=-sk(j)
+ j=j+1
+ if(j .lt. k) go to 630
+ 640 k1=kk
+ k2=kk+kspnn
+ aa=a(kk)
+ bb=b(kk)
+ ak=aa
+ bk=bb
+ j=1
+ k1=k1+kspan
+ 650 k2=k2-kspan
+ j=j+1
+ at(j)=a(k1)+a(k2)
+ ak=at(j)+ak
+ bt(j)=b(k1)+b(k2)
+ bk=bt(j)+bk
+ j=j+1
+ at(j)=a(k1)-a(k2)
+ bt(j)=b(k1)-b(k2)
+ k1=k1+kspan
+ if(k1 .lt. k2) go to 650
+ a(kk)=ak
+ b(kk)=bk
+ k1=kk
+ k2=kk+kspnn
+ j=1
+ 660 k1=k1+kspan
+ k2=k2-kspan
+ jj=j
+ ak=aa
+ bk=bb
+ aj=0.0
+ bj=0.0
+ k=1
+ 670 k=k+1
+ ak=at(k)*ck(jj)+ak
+ bk=bt(k)*ck(jj)+bk
+ k=k+1
+ aj=at(k)*sk(jj)+aj
+ bj=bt(k)*sk(jj)+bj
+ jj=jj+j
+ if(jj .gt. jf) jj=jj-jf
+ if(k .lt. jf) go to 670
+ k=jf-j
+ a(k1)=ak-bj
+ b(k1)=bk+aj
+ a(k2)=ak+bj
+ b(k2)=bk-aj
+ j=j+1
+ if(j .lt. k) go to 660
+ kk=kk+kspnn
+ if(kk .le. nn) go to 640
+ kk=kk-nn
+ if(kk .le. kspan) go to 640
+! multiply by rotation factor (except for factors of 2 and 4)
+ 700 if(i .eq. m) go to 800
+ kk=jc+1
+ 710 c2=1.0-cd
+ s1=sd
+ 720 c1=c2
+ s2=s1
+ kk=kk+kspan
+ 730 ak=a(kk)
+ a(kk)=c2*ak-s2*b(kk)
+ b(kk)=s2*ak+c2*b(kk)
+ kk=kk+kspnn
+ if(kk .le. nt) go to 730
+ ak=s1*s2
+ s2=s1*c2+c1*s2
+ c2=c1*c2-ak
+ kk=kk-nt+kspan
+ if(kk .le. kspnn) go to 730
+ c2=c1-(cd*c1+sd*s1)
+ s1=s1+(sd*c1-cd*s1)
+ c1=2.0-(c2**2+s1**2)
+ s1=c1*s1
+ c2=c1*c2
+ kk=kk-kspnn+jc
+ if(kk .le. kspan) go to 720
+ kk=kk-kspan+jc+inc
+ if(kk .le. jc+jc) go to 710
+ go to 100
+! permute the results to normal order---done in two stages
+! permutation for square factors of n
+ 800 np(1)=ks
+ if(kt .eq. 0) go to 890
+ k=kt+kt+1
+ if(m .lt. k) k=k-1
+ j=1
+ np(k+1)=jc
+ 810 np(j+1)=np(j)/nfac(j)
+ np(k)=np(k+1)*nfac(j)
+ j=j+1
+ k=k-1
+ if(j .lt. k) go to 810
+ k3=np(k+1)
+ kspan=np(2)
+ kk=jc+1
+ k2=kspan+1
+ j=1
+ if(n .ne. ntot) go to 850
+! permutation for single-variate transform (optional code)
+ 820 ak=a(kk)
+ a(kk)=a(k2)
+ a(k2)=ak
+ bk=b(kk)
+ b(kk)=b(k2)
+ b(k2)=bk
+ kk=kk+inc
+ k2=kspan+k2
+ if(k2 .lt. ks) go to 820
+ 830 k2=k2-np(j)
+ j=j+1
+ k2=np(j+1)+k2
+ if(k2 .gt. np(j)) go to 830
+ j=1
+ 840 if(kk .lt. k2) go to 820
+ kk=kk+inc
+ k2=kspan+k2
+ if(k2 .lt. ks) go to 840
+ if(kk .lt. ks) go to 830
+ jc=k3
+ go to 890
+! permutation for multivariate transform
+ 850 k=kk+jc
+ 860 ak=a(kk)
+ a(kk)=a(k2)
+ a(k2)=ak
+ bk=b(kk)
+ b(kk)=b(k2)
+ b(k2)=bk
+ kk=kk+inc
+ k2=k2+inc
+ if(kk .lt. k) go to 860
+ kk=kk+ks-jc
+ k2=k2+ks-jc
+ if(kk .lt. nt) go to 850
+ k2=k2-nt+kspan
+ kk=kk-nt+jc
+ if(k2 .lt. ks) go to 850
+ 870 k2=k2-np(j)
+ j=j+1
+ k2=np(j+1)+k2
+ if(k2 .gt. np(j)) go to 870
+ j=1
+ 880 if(kk .lt. k2) go to 850
+ kk=kk+jc
+ k2=kspan+k2
+ if(k2 .lt. ks) go to 880
+ if(kk .lt. ks) go to 870
+ jc=k3
+ 890 if(2*kt+1 .ge. m) return
+ kspnn=np(kt+1)
+! permutation for square-free factors of n
+ j=m-kt
+ nfac(j+1)=1
+ 900 nfac(j)=nfac(j)*nfac(j+1)
+ j=j-1
+ if(j .ne. kt) go to 900
+ kt=kt+1
+ nn=nfac(kt)-1
+ if(nn .gt. maxp) go to 998
+ jj=0
+ j=0
+ go to 906
+ 902 jj=jj-k2
+ k2=kk
+ k=k+1
+ kk=nfac(k)
+ 904 jj=kk+jj
+ if(jj .ge. k2) go to 902
+ np(j)=jj
+ 906 k2=nfac(kt)
+ k=kt+1
+ kk=nfac(k)
+ j=j+1
+ if(j .le. nn) go to 904
+! determine the permutation cycles of length greater than 1
+ j=0
+ go to 914
+ 910 k=kk
+ kk=np(k)
+ np(k)=-kk
+ if(kk .ne. j) go to 910
+ k3=kk
+ 914 j=j+1
+ kk=np(j)
+ if(kk .lt. 0) go to 914
+ if(kk .ne. j) go to 910
+ np(j)=-j
+ if(j .ne. nn) go to 914
+ maxf=inc*maxf
+! reorder a and b, following the permutation cycles
+ go to 950
+ 924 j=j-1
+ if(np(j) .lt. 0) go to 924
+ jj=jc
+ 926 kspan=jj
+ if(jj .gt. maxf) kspan=maxf
+ jj=jj-kspan
+ k=np(j)
+ kk=jc*k+ii+jj
+ k1=kk+kspan
+ k2=0
+ 928 k2=k2+1
+ at(k2)=a(k1)
+ bt(k2)=b(k1)
+ k1=k1-inc
+ if(k1 .ne. kk) go to 928
+ 932 k1=kk+kspan
+ k2=k1-jc*(k+np(k))
+ k=-np(k)
+ 936 a(k1)=a(k2)
+ b(k1)=b(k2)
+ k1=k1-inc
+ k2=k2-inc
+ if(k1 .ne. kk) go to 936
+ kk=k2
+ if(k .ne. j) go to 932
+ k1=kk+kspan
+ k2=0
+ 940 k2=k2+1
+ a(k1)=at(k2)
+ b(k1)=bt(k2)
+ k1=k1-inc
+ if(k1 .ne. kk) go to 940
+ if(jj .ne. 0) go to 926
+ if(j .ne. 1) go to 924
+ 950 j=k3+1
+ nt=nt-kspnn
+ ii=nt-inc+1
+ if(nt .ge. 0) go to 924
+ return
+! error finish, insufficient array storage
+ 998 isn=0
+ print*,'array bounds exceeded within subroutine fft'
+
+ stop
+ end
Property changes on: DART/branches/development/models/sqg/fft.f90
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/models/sqg/model_description
===================================================================
--- DART/branches/development/models/sqg/model_description (rev 0)
+++ DART/branches/development/models/sqg/model_description 2012-05-31 22:47:55 UTC (rev 5743)
@@ -0,0 +1,42 @@
+
+This is a uniform PV two-surface QG+1 spectral model contributed
+by Rahul Majahan.
+
+The underlying model is described in:
+Hakim, Gregory J., 2000: Role of Nonmodal Growth and Nonlinearity in
+Cyclogenesis Initial-Value Problems. J. Atmos. Sci., 57, 2951–2967.
+
+
+NOTE:
+
+This model counts on a non-standard interface to the obs_model_mod.f90
+code in DART/obs_model. To compile and run this model, make the
+following (non-backwards-compatible) two changes:
+
+
+1:
+
+ !------------- Block for subroutine callable adv_1step interface -----------
+ if(async == 0) then
+
+ do while(ens_handle%time(i) < target_time)
+ call adv_1step(ens_handle%vars(:, i), ens_handle%time(i))
+ ens_handle%time(i) = ens_handle%time(i) + time_step
+ end do
+
+! start of new section
+ else if(async == 1) then
+
+ call adv_1step(ens_handle%vars(:, i), ens_handle%time(i), target_time)
+ ens_handle%time(i) = target_time
+! end of new section
+
+ !-------------- End single subroutine callable adv_1step interface ---------
+
+and 2:
+
+! Following is for async options that use shell to advance model
+SHELL_ADVANCE_METHODS: if(async > 1) then
+
+! changed from: if(async /= 0) then
+
Added: DART/branches/development/models/sqg/model_mod.f90
===================================================================
--- DART/branches/development/models/sqg/model_mod.f90 (rev 0)
+++ DART/branches/development/models/sqg/model_mod.f90 2012-05-31 22:47:55 UTC (rev 5743)
@@ -0,0 +1,1057 @@
+
+! DART software - Copyright 2004 - 2011 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+module model_mod
+
+!========================================================================
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+!========================================================================
+
+!========================================================================
+! Assimilation interface for:
+! Uniform-PV two-surface QG+1 model in spectral form (Hakim 2000)
+!========================================================================
+
+!================= m o d u l e i n f o r m a t i o n ==================
+
+use types_mod, only : r8, MISSING_R8
+use time_manager_mod, only : time_type, set_time, get_time, &
+ operator(<), operator(+)
+use location_mod, only : location_type, get_close_maxdist_init, &
+ get_close_obs_init, get_close_obs, set_location, &
+ get_location, set_location_missing
+use utilities_mod, only : register_module, error_handler, nc_check, &
+ E_ERR, E_MSG, logfileunit, get_unit, close_file, &
+ dump_unit_attributes, &
+ nmlfileunit, do_output, do_nml_file, do_nml_term, &
+ find_namelist_in_file, check_namelist_read
+use obs_kind_mod, only : KIND_POTENTIAL_TEMPERATURE
+use spectral_mod
+use sqg_mod, only : diffusion, init, init_jet, terrain, invert, advect, &
+ tadv, xy_to_sp, sp_to_xy, d_setup, ft_2d, norm
+
+!========================================================================
+
+implicit none
+
+private
+
+public :: get_model_size, &
+ adv_1step, &
+ get_state_meta_data, &
+ model_interpolate, &
+ get_model_time_step, &
+ end_model, &
+ static_init_model, &
+ init_time, &
+ init_conditions, &
+ nc_write_model_atts, &
+ nc_write_model_vars, &
+ pert_model_state, &
+ get_close_maxdist_init, &
+ get_close_obs_init, &
+ get_close_obs, &
+ ens_mean_for_model
+
+! public subroutines / functions
+public :: sqg_to_dart, &
+ dart_to_sqg, &
+ get_model_static_data
+
+! public types
+public :: model_static
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+!---------------------------------------------------------------
+! Namelist with default values
+
+logical :: output_state_vector = .false.
+real(r8) :: channel_center = 45.0_r8
+real(r8) :: channel_width = 40.0_r8
+logical :: debug = .false.
+namelist /model_nml/ output_state_vector, &
+ channel_center, channel_width, &
+ debug
+
+real, dimension(mmax,nmax), parameter :: Rblank = 0.0
+complex, dimension(mmax,nmax), parameter :: Cblank = 0.0
+
+integer, parameter :: model_size = 2 * (2*kmax) * (2*lmax)
+
+! Define the location of the state variables in module storage
+type(location_type), allocatable :: state_loc(:)
+
+type(time_type) :: time_step
+character(len=129) :: errstring
+
+type model_static
+ logical :: top, bot
+ real :: dco, lam
+ complex, allocatable, dimension(:,:) :: thbB, thbT
+ real, allocatable, dimension(:,:) :: thbyB, thbyT, ulinB, ulinT
+ real, allocatable, dimension(:,:) :: hx, hy, hu, hv
+ real(r8), allocatable, dimension(:) :: lons, lats, levs
+end type model_static
+
+type(model_static) :: sqg_static
+
+!========================================================================
+
+contains
+
+!========================================================================
+
+subroutine static_init_model()
+!-----------------------------
+! Initializes class data for surface quasigeostrophy model
+! (all the stuff that needs to be done once.)
+
+integer :: iunit, io
+
+integer :: i, lev_index, lat_index, lon_index
+real(r8) :: lon, lat, lev
+
+! Print module information to log file and stdout.
+call register_module(source, revision, revdate)
+
+! This is where you would read a namelist, for example.
+call find_namelist_in_file("input.nml", "model_nml", iunit)
+read(iunit, nml = model_nml, iostat = io)
+call check_namelist_read(iunit, io, "model_nml")
+
+! Record the namelist values used for the run ...
+if (do_nml_file()) write(nmlfileunit, nml=model_nml)
+if (do_nml_term()) write( * , nml=model_nml)
+
+! evenly spaced longitudes, latitudes and two levels
+allocate( sqg_static%lons(2*kmax) )
+allocate( sqg_static%lats(2*lmax) )
+allocate( sqg_static%levs(2 ) )
+do i = 1, 2*kmax
+ sqg_static%lons(i) = 360.0*(i-1.0)/(2*kmax)
+enddo
+! channel is centered about "channel_center" and has a width of "channel_width"
+do i = 1, 2*lmax
+ sqg_static%lats(i) = channel_center + channel_width * (real(i,r8)/(2*lmax) - 0.5)
+enddo
+do i = 1, 2
+ sqg_static%levs(i) = i
+enddo
+
+! allocate storage for locations
+allocate( state_loc(model_size) )
+do i = 1, model_size
+ lev_index = (i-1) / ((2*lmax)*(2*kmax)) + 1
+ lat_index = ((i-1) - ((lev_index-1)*(2*lmax)*(2*kmax))) / (2*kmax) + 1
+ lon_index = (i-1) - ((lev_index-1)*(2*lmax)*(2*kmax)) - ((lat_index-1)*(2*kmax)) + 1
+
+ lev = sqg_static%levs(lev_index)
+ lat = sqg_static%lats(lat_index)
+ lon = sqg_static%lons(lon_index)
+
+ ! With the location module ... you specify that the
+ ! vertical coordinate is a 'level' by 'which_vert' == 1
+ state_loc(i) = set_location(lon, lat, lev, 1)
+enddo
+
+! get dimensional model time-step,
+! follow non-dimensionalization conversion for sQG (Mahajan & Hakim 2011, Table 1)
+! dTs = (Ls/Us) * dt ; Ls = 1000 km = 1e6 m ; Us = 30 m/s
+time_step = set_time(int((1e6/30)*dt), 0)
+
+! initialize derivative operators in the model
+call d_setup()
+
+! initialize diffusion coefficient
+call diffusion(sqg_static%dco)
+
+! advection flags
+if ( (model .eq. 0) .or. (model .eq. 1) ) then ! 2D or 2sQG
+ sqg_static%top = .true.; sqg_static%bot = .true.
+elseif ( model .eq. 2 ) then ! tropo sQG
+ sqg_static%top = .true.; sqg_static%bot = .false.
+elseif ( (model .eq. 3) .or. (model .eq. 0) ) then ! surface sQG
+ sqg_static%top = .false.; sqg_static%bot = .true.
+elseif ( model .eq. 4 ) then ! tropo HsQG
+ sqg_static%top = .true.; sqg_static%bot = .false.
+endif
+!if (maxval(abs(thxyT)) .lt. 1.e-5) sqg_static%top = .false.
+!if (maxval(abs(thxyB)) .lt. 1.e-5) sqg_static%bot = .false.
+
+! initialize Hoskins-West jet
+allocate(sqg_static%thbB(2*kmax,2*lmax)) ; allocate(sqg_static%thbT(2*kmax,2*lmax))
+allocate(sqg_static%thbyB(mmax,nmax)) ; allocate(sqg_static%thbyT(mmax,nmax))
+allocate(sqg_static%ulinB(mmax,nmax)) ; allocate(sqg_static%ulinT(mmax,nmax))
+if ( hw ) then
+
+ call init_jet(sqg_static%thbB, sqg_static%thbT, &
+ sqg_static%thbyB, sqg_static%thbyT, &
+ sqg_static%ulinB, sqg_static%ulinT, &
+ sqg_static%lam)
+
+ ! barotropic wind (Ross = 0!)
+ sqg_static%ulinB = sqg_static%ulinB + 0.0
+ sqg_static%ulinT = sqg_static%ulinT - 0.0 * H * sqg_static%lam
+
+else
+
+ sqg_static%thbB = 0.0 ; sqg_static%thbT = 0.0
+ sqg_static%thbyB = 0.0 ; sqg_static%thbyT = 0.0
+ sqg_static%ulinB = 0.0 ; sqg_static%ulinT = 0.0
+ sqg_static%lam = 0.0
+
+endif
+
+! initialize terrain
+allocate(sqg_static%hx(mmax,nmax)) ; allocate(sqg_static%hy(mmax,nmax))
+allocate(sqg_static%hu(mmax,nmax)) ; allocate(sqg_static%hv(mmax,nmax))
+if ( (iterr) .and. (Ross .eq. 0.0) ) then
+ call terrain(sqg_static%hx, sqg_static%hy, &
+ sqg_static%hu, sqg_static%hv)
+else
+ sqg_static%hx = 0.0 ; sqg_static%hy = 0.0
+ sqg_static%hu = 0.0 ; sqg_static%hv = 0.0
+endif
+
+end subroutine static_init_model
+
+!========================================================================
+
+subroutine init_conditions(x)
+!----------------------------
+! read initial conditions from a file,
+! must only contain perturbation theta
+
+real(r8), intent(out) :: x(:)
+
+integer :: j
+complex, allocatable, dimension(:,:) :: thspB, thspT
+real, allocatable, dimension(:,:) :: thxyB, thxyT
+real, allocatable, dimension(:,:,:) :: theta
+
+! initialize the theta fields:
+allocate(thxyB(2*kmax,2*lmax)) ; allocate(thxyT(2*kmax,2*lmax))
+call init('sqgInput.nc',thxyB,thxyT)
+
+! first add base-state jet:
+if ( hw ) then
+
+ allocate(thspB(2*kmax,2*lmax)) ; allocate(thspT(2*kmax,2*lmax))
+
+ ! map into spectral space at the same resolution:
+ call xy_to_sp(cmplx(thxyB,0.0),thspB,2*kmax,2*lmax,kmax,lmax)
+ call xy_to_sp(cmplx(thxyT,0.0),thspT,2*kmax,2*lmax,kmax,lmax)
+
+ thspB = thspB + sqg_static%thbB
+ thspT = thspT + sqg_static%thbT
+
+ ! map into grid-point space space at the same resolution:
+ call sp_to_xy(thspB,thxyB,kmax,lmax,2*kmax,2*lmax)
+ call sp_to_xy(thspT,thxyT,kmax,lmax,2*kmax,2*lmax)
+
+ deallocate(thspB) ; deallocate(thspT)
+
+endif
+
+! second add linear shear:
+allocate(theta(2*kmax,2*lmax,2))
+do j = 1, 2*lmax
+ theta(:,j,1) = thxyB(:,j) - sqg_static%lam * real(j-1) * YL/real(2*lmax)
+ theta(:,j,2) = thxyT(:,j) - sqg_static%lam * real(j-1) * YL/real(2*lmax)
+enddo
+
+! wrap model variables into DART state vector:
+call sqg_to_dart(theta,x)
+
+deallocate(thxyB) ; deallocate(thxyT)
+deallocate(theta)
+
+end subroutine init_conditions
+
+!========================================================================
+
+subroutine adv_1step(x, current_time, target_time)
+!-------------------------------------------------
+! Does a time advance for sQG model with state vector as
+! input and output.
+! This interface advances from current_time to target_time
+! and should be used with async = -2
+! see modifications in obs_model_mod.f90
+
+real(r8), intent(inout) :: x(:)
+type(time_type), intent(in) :: current_time
+type(time_type), intent(in), optional :: target_time
+
+real(r8) :: cxB,cyB,cxT,cyT
+integer :: j
+integer :: itime, cdays, cseconds, tdays, tseconds
+type(time_type) :: ctime
+
+real, allocatable, dimension(:,:,:) :: thxy
+real, allocatable, dimension(:,:) :: thxyB, thxyT
+real, allocatable, dimension(:,:) :: laplacian
+real, allocatable, dimension(:,:) :: thxB, thxT, thyB, thyT
+real, allocatable, dimension(:,:) :: uB, uT, vB, vT
+real, allocatable, dimension(:,:) :: tthB, tthT
+
+complex, allocatable, dimension(:,:) :: thspB, thspT, thspB1, thspT1, thspB2, thspT2
+complex, allocatable, dimension(:,:) :: sB, sBold
+complex, allocatable, dimension(:,:) :: tthspB, tthspT
+
+logical :: first
+
+! allocate necessary space up-front:
+allocate(thxy(2*kmax,2*lmax,2))
+allocate(thxyB( 2*kmax,2*lmax)) ; allocate(thxyT( 2*kmax,2*lmax))
+
+allocate(laplacian(mmax,nmax))
+allocate(thxB(mmax,nmax)) ; allocate(thxT(mmax,nmax))
+allocate(thyB(mmax,nmax)) ; allocate(thyT(mmax,nmax))
+allocate(uB( mmax,nmax)) ; allocate(uT( mmax,nmax))
+allocate(vB( mmax,nmax)) ; allocate(vT( mmax,nmax))
+allocate(tthB(mmax,nmax)) ; allocate(tthT(mmax,nmax))
+
+allocate(thspB( 2*kmax,2*lmax)) ; allocate(thspT( 2*kmax,2*lmax))
+allocate(thspB1(2*kmax,2*lmax)) ; allocate(thspT1(2*kmax,2*lmax))
+allocate(thspB2(2*kmax,2*lmax)) ; allocate(thspT2(2*kmax,2*lmax))
+
+allocate(sB(2*kmax,2*lmax)) ; allocate(sBold(2*kmax,2*lmax))
+
+allocate(tthspB(mmax,nmax)) ; allocate(tthspT(mmax,nmax))
+
+! unwrap DART state vector into model variables:
+call dart_to_sqg(x, thxy)
+
+! first remove the linear shear:
+do j = 1, 2*lmax
+ thxyB(:,j) = thxy(:,j,1) + sqg_static%lam * real(j-1) * YL/real(2*lmax)
+ thxyT(:,j) = thxy(:,j,2) + sqg_static%lam * real(j-1) * YL/real(2*lmax)
+enddo
+
+! map into spectral space at the same resolution:
+call xy_to_sp(cmplx(thxyB,0.0), thspB, 2*kmax,2*lmax, kmax, lmax)
+call xy_to_sp(cmplx(thxyT,0.0), thspT, 2*kmax,2*lmax, kmax, lmax)
+
+! second remove the HW base-state:
+if ( hw ) then
+ thspB = thspB - sqg_static%thbB
+ thspT = thspT - sqg_static%thbT
+endif
+
+! initialize certain variables for first time-step
+first = .true.
+sB = 0.0
+thspB1 = 0.0; thspB2 = 0.0
+thspT1 = 0.0; thspT2 = 0.0
+
+! advance from current_time to target_time:
+itime = 1
+ctime = current_time
+do while ( ctime < target_time )
+
+ ! save old stream-function for Ekman calculation:
+ sBold = sB
+
+ ! invert theta for stream-function; compute gradients for advection:
+ call invert(thspB, thspT, thxB, thxT, thyB, thyT, vB, vT, uB, uT, &
+ sqg_static%thbB, sqg_static%thbT, sqg_static%thbyB, sqg_static%thbyT, &
+ sqg_static%ulinB, sqg_static%ulinT, &
+ first, sqg_static%bot, sqg_static%top, sqg_static%lam, &
+ sB, sBold, laplacian)
+
+ ! option to compute potential enstrophy norm and growth rate:
+ if ( inorm ) call norm(thspB, thspT, itime)
+
+ ! spectral advection:
+ if ( sqg_static%bot ) &
+ call advect(uB, vB, thxB, thyB, &
+ sqg_static%thbyB, sqg_static%hx, sqg_static%hy, sqg_static%ulinB, &
+ tthB, sqg_static%lam, laplacian)
+ if ( sqg_static%top ) &
+ call advect(uT + sqg_static%hu, vT + sqg_static%hv, thxT, thyT, &
+ sqg_static%thbyT, Rblank, Rblank, sqg_static%ulinT, &
+ tthT, sqg_static%lam, Rblank)
+
+ ! compute Courant numbers and print to stdout if debugging:
+ if ( mod((itime-1),10) .eq. 0 ) then
+ cxB = maxval(abs(uB + sqg_static%ulinB)) * dt / (XL/real(2*kmax))
+ cyB = maxval(abs(vB)) * dt / (YL/real(2*lmax))
+ cxT = maxval(abs(uT + sqg_static%ulinT)) * dt / (XL/real(2*kmax))
+ cyT = maxval(abs(vT)) * dt / (YL/real(2*lmax))
+ if ( debug ) write(*,'(A23,F10.3,4F8.3))') 'time,cxB,cyB,cxT,cyT = ', &
+ real(itime-1)*dt,cxB,cyB,cxT,cyT
+ endif
+
+ ! FFT back to spectral space:
+ if ( sqg_static%bot ) then
+ tthspB = cmplx(tthB,0.0)
+ call ft_2d(tthspB,mmax,nmax,-1)
+ endif
+ if ( sqg_static%top ) then
+ tthspT = cmplx(tthT,0.0)
+ call ft_2d(tthspT,mmax,nmax,-1)
+ endif
+
+ ! advance one time-step with explicit (hyper-) diffusion:
+ if ( sqg_static%bot ) call tadv(thspB, tthspB, thspB1, thspB2, sqg_static%dco, first)
+ if ( sqg_static%top ) call tadv(thspT, tthspT, thspT1, thspT2, sqg_static%dco, first)
+
+ ! zero out k=K, l=L modes on the bottom boundary:
+ thspB(kmax,:) = 0.0 ; thspB(lmax,:) = 0.0
+
+ itime = itime + 1
+ ctime = ctime + time_step
+
+ first = .false.
+
+enddo
+
+! first add the HW base-state:
+if ( hw ) then
+ thspB = thspB + sqg_static%thbB
+ thspT = thspT + sqg_static%thbT
+endif
+
+! map into grid-point space at the same resolution:
+call sp_to_xy(thspB, thxyB, kmax, lmax, 2*kmax, 2*lmax)
+call sp_to_xy(thspT, thxyT, kmax, lmax, 2*kmax, 2*lmax)
+
+! second add the linear shear:
+do j = 1, 2*lmax
+ thxy(:,j,1) = thxyB(:,j) - sqg_static%lam * real(j-1) * YL/real(2*lmax)
+ thxy(:,j,2) = thxyT(:,j) - sqg_static%lam * real(j-1) * YL/real(2*lmax)
+enddo
+
+! wrap model variables into DART state vector:
+call sqg_to_dart(thxy, x)
+
+! deallocate space allocated up-front:
+deallocate(thxy)
+deallocate(thxyB ) ; deallocate(thxyT )
+
+deallocate(laplacian)
+deallocate(thxB) ; deallocate(thxT)
+deallocate(thyB) ; deallocate(thyT)
+deallocate(uB) ; deallocate(uT)
+deallocate(vB) ; deallocate(vT)
+deallocate(tthB) ; deallocate(tthT)
+
+deallocate(thspB ) ; deallocate(thspT )
+deallocate(thspB1) ; deallocate(thspT1)
+deallocate(thspB2) ; deallocate(thspT2)
+
+deallocate(sB) ; deallocate(sBold)
+
+deallocate(tthspB) ; deallocate(tthspT)
+
+end subroutine adv_1step
+
+!========================================================================
+
+function get_model_size()
+!------------------------
+! Returns the size of the model as an integer.
+
+integer :: get_model_size
+
+get_model_size = model_size
+
+end function get_model_size
+
+!========================================================================
+
+subroutine init_time(time)
+!-------------------------
+! For now, returns the initialization time as 0
+
+type(time_type), intent(out) :: time
+
+! for now, just set to 0
+print*,'init_time'
+time = set_time(0,0)
+
+end subroutine init_time
+
+!========================================================================
+
+subroutine model_interpolate(x, location, itype, obs_val, istatus)
+!-----------------------------------------------------------------
+! Interpolates from state vector x to the location and returns obs_val.
+! istatus = 0 suggests interpolation went ok, 1 means something went wrong.
+! Argument itype specifies the type of field (for instance potential temperature in this model)
+
+real(r8), intent(in) :: x(:)
+type(location_type), intent(in) :: location
+integer, intent(in) :: itype
+real(r8), intent(out) :: obs_val
+integer, intent(out) :: istatus
+
+real(r8) :: lon, lat, lev, lon_lat_lev(3)
+real(r8) :: bot_lon, top_lon, delta_lon, bot_lat, top_lat, delta_lat
+real(r8) :: temp_lon, lon_fract, lat_fract, val(2,2), a(2)
+integer :: lon_below, lon_above, lat_below, lat_above, level, i
+
+! Assume all interpolations okay for now
+istatus = 0
+
+lon_lat_lev = get_location(location)
+lon = lon_lat_lev(1)
+lat = lon_lat_lev(2)
+lev = lon_lat_lev(3); level = int(lev)
+
+! Level is obvious for now, but make sure that it is within valid range
+if ( (level < 1) .or. (level > 2) ) then
+ istatus = 1
+ obs_val = MISSING_R8
+ write(errstring,*)'level ',level,' must be 1 or 2'
+ call error_handler(E_ERR,'model_mod:model_interpolate', errstring, source, revision, revdate)
+endif
+
+! Get globally defined lat / lon grid specs
+bot_lon = sqg_static%lons(1)
+top_lon = sqg_static%lons(2*kmax)
+bot_lat = sqg_static%lats(1)
+top_lat = sqg_static%lats(2*lmax)
+delta_lon = sqg_static%lons(2) - sqg_static%lons(1)
+delta_lat = sqg_static%lats(2) - sqg_static%lats(1)
+
+! Compute bracketing lon indices
+if (lon >= bot_lon .and. lon <= top_lon ) then
+ lon_below = int((lon - bot_lon) / delta_lon) + 1
+ lon_above = lon_below + 1
+ lon_fract = (lon - ((lon_below - 1) * delta_lon + bot_lon)) / delta_lon
+else
+ ! At wraparound point
+ lon_below = 2*kmax
+ lon_above = 1
+ if(lon < bot_lon) then
+ temp_lon = lon + 360.0_r8
+ else
+ temp_lon = lon
+ endif
+ lon_fract = (temp_lon - top_lon) / delta_lon
+endif
+
+! Next, compute bracketing lat indices
+if (lat >= bot_lat .and. lat <= top_lat) then
+ lat_below = int((lat - bot_lat) / delta_lat) + 1
+ lat_above = lat_below + 1
+ lat_fract = (lat - ((lat_below - 1) * delta_lat + bot_lat)) / delta_lat
+else
+ ! Outside the channel, return with an error status
+ istatus = 1
+ obs_val = MISSING_R8
+ return
+endif
+
+! Get values at the 4 surrounding points
+val(1,1) = get_val(x, lon_below, lat_below, level)
+val(1,2) = get_val(x, lon_below, lat_above, level)
+val(2,1) = get_val(x, lon_above, lat_below, level)
+val(2,2) = get_val(x, lon_above, lat_above, level)
+
+! Do the weighted average for interpolation
+if ( debug ) write(*,*) 'fracts ', lon_fract, lat_fract
+do i = 1, 2
+ a(i) = lon_fract * val(2, i) + (1.0 - lon_fract) * val(1, i)
+end do
+
+obs_val = lat_fract * a(2) + (1.0 - lat_fract) * a(1)
+
+end subroutine model_interpolate
+
+!========================================================================
+
+function get_val(x, lon_index, lat_index, level)
+!-----------------------------------------------
+! Given the state vector and location, returns the value at that location
+
+real(r8), intent(in) :: x(:)
+integer, intent(in) :: lon_index, lat_index, level
+real(r8) :: get_val
+
+integer :: indx
+
+indx = (level-1)*(2*lmax)*(2*kmax) + (lat_index-1)*(2*kmax) + lon_index
+
+!! should not be possible now; but this error check can be commented back in.
+!! (it is out for performance reasons, but if you get any strange values, this
+!! is a good first check to re-enable.)
+!if (indx < 1 .or. indx > size(x)) then
+! write(errstring,*)'index ',indx,' not between 1 and ', size(x), ' (should not be possible)'
+! call error_handler(E_ERR,'model_mod:get_val', errstring, source, revision, revdate)
+!endif
+
+get_val = x(indx)
+
+end function get_val
+
+!========================================================================
+
+function get_model_time_step()
+!-----------------------------
+! Returns the time step of the model
+
+type(time_type) :: get_model_time_step
+
+get_model_time_step = time_step
+
+end function get_model_time_step
+
+!========================================================================
+
+subroutine get_state_meta_data(index_in, location, var_type)
+!-----------------------------------------------------------
+! Given an integer index into the state vector structure, returns the
+! associated location.
+
+integer, intent(in) :: index_in
+type(location_type), intent(out) :: location
+integer, intent(out), optional :: var_type
+
+integer :: lev_index, lat_index, lon_index
+real(r8) :: lon, lat, lev
+integer :: location_index
+
+! avoid out-of-range queries
+if ( index_in > model_size ) then
+ write(errstring,*)'index_in ',index_in,' must be between 1 and ', model_size
+ call error_handler(E_ERR,'model_mod:get_state_meta_data', errstring, source, revision, revdate)
+endif
+
+lev_index = (index_in-1) / ((2*lmax)*(2*kmax)) + 1
+lat_index = ((index_in-1) - ((lev_index-1)*(2*lmax)*(2*kmax))) / (2*kmax) + 1
+lon_index = (index_in-1) - ((lev_index-1)*(2*lmax)*(2*kmax)) - ((lat_index-1)*(2*kmax)) + 1
+
+lev = sqg_static%levs(lev_index)
+lat = sqg_static%lats(lat_index)
+lon = sqg_static%lons(lon_index)
+
+! With the threed_sphere location module ... you specify that the
+! vertical coordinate is a 'level' by 'which_vert' == 1
+location = set_location(lon, lat, lev, 1)
+
+! Alternately, use state_loc that is defined in static_init_model()
+!location = state_loc(index_in)
+
+if (present(var_type)) var_type = KIND_POTENTIAL_TEMPERATURE
+
+end subroutine get_state_meta_data
+
+!========================================================================
+
+subroutine end_model()
+!---------------------
+! Does any shutdown and clean-up needed for model.
+
+! good practice ... deallocate stuff from static_init_model
+deallocate(sqg_static%lons) ; deallocate(sqg_static%lats) ; deallocate(sqg_static%levs)
+
+deallocate(sqg_static%thbB ) ; deallocate(sqg_static%thbT )
+deallocate(sqg_static%thbyB) ; deallocate(sqg_static%thbyT)
+deallocate(sqg_static%ulinB) ; deallocate(sqg_static%ulinT)
+
+deallocate(sqg_static%hx) ; deallocate(sqg_static%hy)
+deallocate(sqg_static%hu) ; deallocate(sqg_static%hv)
+
+deallocate(state_loc)
+
+end subroutine end_model
+
+!========================================================================
+
+function nc_write_model_atts( ncFileID ) result (ierr)
+!-----------------------------------------------------
+! Writes the model-specific attributes to a netCDF file.
+
+use typeSizes
+use netcdf
+
+integer, intent(in) :: ncFileID ! netCDF file identifier
+integer :: ierr ! return value of function
+
+integer :: nDimensions, nVariables, nAttributes, unlimitedDimID
+
+integer :: latDimID, lonDimID, levDimID, StateVarDimID, MemberDimID, TimeDimID
+integer :: latVarID, lonVarID, levVarID, StateVarVarID, StateVarID
+
+character(len=8) :: crdate ! needed by F90 DATE_AND_TIME intrinsic
+character(len=10) :: crtime ! needed by F90 DATE_AND_TIME intrinsic
+character(len=5) :: crzone ! needed by F90 DATE_AND_TIME intrinsic
+integer, dimension(8) :: values ! needed by F90 DATE_AND_TIME intrinsic
+character(len=NF90_MAX_NAME) :: str1
+character(len=NF90_MAX_NAME) :: filename
+
+integer :: i
+
+ierr = -1 ! assume things go poorly
+
+!--------------------------------------------------------------------
+! we only have a netcdf handle here so we do not know the filename
+! or the fortran unit number. but construct a string with at least
+! the netcdf handle, so in case of error we can trace back to see
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list