program testlats C C Sample program to verify LATS is working properly and C to demonstrate a typical AMIP II application C C Mike Fiorino, PCMDI (lats@pcmdi.llnl.gov) C version 1.0 C 10 March, 1997 C C******* WARNING C******* machine / compiler dependent include line C******* C C******* there is no standard but this worked on all plaformts. C******* NB - this is a relative path, and needs to be adjusted C******* depending on where you installed LATS C C include "../include/lats.inc" C C define the grid: C C 72 points in long (no wrap) C 46 points in lat (pole points) C 3 levels C 2 variables C parameter(ni=72,nj=46,nk=3,nv=2) character*20 center character*20 model character*9 var dimension var(nv),id_var(nv) double precision rlon(ni),rlat(nj),plev(nk),slev real ta(ni,nj,nk),psl(ni,nj) C C define the production center and the process (model) C data center/'PCMDI'/ data model/'lats'/ C C set the variable names (AMIP II convention) C C psl - sea level pressure (Pa) C ta - air temperature (degK) C data var/'psl','ta'/ C C define the pressure levels (hPa) C data plev/850.0,500.0,200.0/ C C0000000000000000000000000000000000000000000000000 C C STEP #0 -- set the parameter table C C0000000000000000000000000000000000000000000000000 id_parmtab=latsparmtab("../table/amip2.lats.table") if(id_parmtab.eq.0) stop 'latsparmtab error' C1111111111111111111111111111111111111111111111111 C C STEP #1 -- create a LATS file defining C C 1) convention (e.g., GRIB with a standard calandar) C 2) time statistic (e.g., monthly) C 3) time increment (1 for 1 month) C 4) who or the center (e.g., PCMDI) C 5) what or the model (e.g., lats) C C1111111111111111111111111111111111111111111111111 nconv=2 do iconv=1,nconv if(iconv.eq.1) latsconv=LATS_GRADS_GRIB if(iconv.eq.2) latsconv=LATS_COARDS if(latsconv.eq.LATS_GRADS_GRIB) then id_fil = latscreate('latsout', $ latsconv, $ LATS_STANDARD, $ LATS_MONTHLY,1,center, $ model,'LATS GRIB test') print*,'LATS grib file id = ',id_fil else if(latsconv.eq.LATS_COARDS) then id_fil = latscreate('latsout', $ LATS_COARDS, $ LATS_STANDARD, $ LATS_MONTHLY,1,center, $ model,'LATS netcdf test') print*,'LATS netcdf file id = ',id_fil else go to 800 endif C222222222222222222222222222222222222222222222222222222 C C STEP #2 -- define the grid: C C 1) lons and lats and # points C 2) type (e.g., linear, gaussian or generic) C C do this only once !!! C C222222222222222222222222222222222222222222222222222222 if(iconv.eq.1) then do i=1,ni rlon(i)=0.0+(i-1)*360.0/ni end do do j=1,nj rlat(j)=-90.0+(j-1)*180.0/(nj-1) end do id_grd=latsgrid("u54",LATS_LINEAR, ni, rlon, nj, $ rlat) if(id_grd.eq.0) stop 'latsgrid error' endif C333333333333333333333333333333333333333333333333333333 C C STEP #3 -- define the vertical dimension: C C 1) multi or sfc C 2) levels C C333333333333333333333333333333333333333333333333333333 id_lev=latsvertdim('pressure', 'plev', nk, plev) if(id_lev.eq.0) stop 'latsvertdim error' C444444444444444444444444444444444444444444444444444444 C C STEP #4 -- define the variables C C 1) variable (from the table, e.g., "psl" for mean sea leve pressure) C 2) data type (typically LATS_FLOAT) C 3) time statistic type (e.g., LATS_AVERAGE) C 4) associated vertical dimension from step #3 C C C444444444444444444444444444444444444444444444444444444 id_var(1)=latsvar(id_fil,var(1), $ LATS_FLOAT,LATS_AVERAGE,id_grd, $ 0, 'sfc variable') if(id_var(1).eq.0) stop 'latsvar(1) error' id_var(2)=latsvar(id_fil,var(2), $ LATS_FLOAT,LATS_AVERAGE,id_grd, $ id_lev, 'ua variable') if(id_var(2).eq.0) stop 'latsvar(2) error' C4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a C C STEP #4a (optional) - set a missing value C uncomment the following line C CCCC ierr=latsmissreal(id_fil,id_var(1),1e20,1e13) C C4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a4a C555555555555555555555555555555555555555555555555555555 C C STEP #5 -- get and write the fields C C 1) define the valid time (beginning of the interval) C 2) read(generate) some data C 3) write it out C C555555555555555555555555555555555555555555555555555555 C C valid time 00UTC 1 Jan 1979 C do imo=1,2 iyr=1979 ida=1 ihr=0 C C create sfc field and write C call read_data(psl,var(1),ni,nj,0,imo) C C I know this is overkill, but if you don't do this the hp version C will core dump! C slev=0.0D0 id_write1=latswrite(id_fil,id_var(1),slev, $ iyr,imo,ida,ihr,psl) if(id_write1.eq.0) stop 'latswrite error - sfc' C C create multi-level field C do k=1,nk call read_data(ta(1,1,k),var(2),ni,nj,k,imo) id_write2=latswrite(id_fil,id_var(2),plev(k), $ iyr,imo,ida,ihr, $ ta(1,1,k)) if(id_write2.eq.0) stop 'latswrite error - ua' end do end do C666666666666666666666666666666666666666666666666666666 C C STEP #6 - close the file C C 1) for the GRIB, this creates the GrADS .ctl file for C futher processing by cdunif, GrADS or VCS C C666666666666666666666666666666666666666666666666666666 id_close = latsclose(id_fil) if(id_close.eq.0) stop 'latsclose error' end do C C normal exit C go to 999 C C error conditions C 800 continue print*,'Error 800, undefined LATS output convention' stop 800 999 continue stop end subroutine read_data(a,name,ni,nj,id,it) C C routine to generate quasi realistic AMIP II fields C dimension a(ni,nj) character name*9 pi=3.1459 C C sfc field: psl = sea level pressure C pscl=10.0 pmin=950.0 p0=1013.0 i0=ni/2 j0=nj/2 rly=j0 rlx=i0*0.5 C C low eastward in time C i0=i0+(it-1)*3 if(id.eq.0) then do i=1,ni do j=1,nj x=(i-i0)*1.0 y=(j-j0)*1.0 r=sqrt( x*x + y*y ) p=p0-(p0-pmin)*exp(-r/pscl) a(i,j)=p*100.0 end do end do return endif C C upper air field temperature C C 850 mb C if(id.eq.1) then t0=250.0 ay=10.0 ax=30.0 ishft=0 endif C C 500 mb C if(id.eq.2) then t0=235.0 ay=20.0 ax=20.0 ishft=-5 endif C C 200 mb C if(id.eq.3) then t0=210.0 ay=30.0 ax=10.0 ishft=-10 endif C C shift pattern eastward in time C ishft=ishft-(it-1)*2 do j=1,nj y=j0-(j-1) do i=1,ni angy=pi*0.5*(y/rly) angyx=angy*2 x=i0-(i-1)+ishft angx=pi*0.5*(x/rlx) a(i,j)=t0 + ay*cos(angy) + ax*sin(angyx)*sin(angx) end do end do return end