Index: preproc/src/preconvert.f90
===================================================================
--- preproc/src/preconvert.f90	(revision 57637935b8c736ecb0a376df5e8f6a044bbd05af)
+++ preproc/src/preconvert.f90	(revision 57637935b8c736ecb0a376df5e8f6a044bbd05af)
@@ -0,0 +1,470 @@
+PROGRAM PRECONVERT
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !                                                                 !
+  ! PROGRAM PRECONVERT - PREPARES INPUT DATA FOR POP MODEL METEOR-  !
+  !                      OLOGICAL PREPROCESSOR                      !
+  !                                                                 !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !                                                                 !
+  ! CALCULATION OF ETAPOINT ON A REGULAR LAMDA/PHI GRID AND WRITING !
+  ! U,V,ETAPOINT,T,PS,Q,SD,MSL,TCC,10U, 10V, 2T,2D,LSP,CP,SSHF,SSR, !
+  ! EWSS,NSSS TO AN OUTPUT FILE (GRIB 1 or 2 FORMAT).               ! 
+  !                                                                 !
+  ! AUTHORS: L. HAIMBERGER, G. WOTAWA, 1994-04                      !
+  !                     adapted: A. BECK                            !
+  !                     2003-05-11                                  !
+  !          L. Haimberger 2006-12    V2.0                          !
+  !                    modified to handle arbitrary regular grids   !
+  !                    and T799 resolution data                     !
+  !          L. Haimberger 2010-03    V4.0                          !
+  !                    modified to grib edition 2 fields            !
+  !                    and T1279 resolution data                    !
+  !                                                                 !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !                                                                 !
+  ! DESCRIPTION OF NEEDED INPUT:                                    !
+  !                                                                 !
+  !       FILE      PARAMETER(S)    DATA REPRESENTATION             !
+  !                                                                 !
+  !       fort.10   U,V             spherical harmonics             !
+  !       fort.11   T               regular lamda phi grid          !
+  !       fort.12   LNSP            spherical harmonics             !
+  !       fort.13   D               spherical harmonics             !
+  !       fort.14   SD,MSL,TCC,10U,                                 !
+  !                 10V,2T,2D       regular lamda phi grid          !
+  !       fort.18   Q               regular lamda phi grid          !
+  !                                                                 !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !                                                                 !
+  ! DESCRIPTION OF OUTPUT:                                          !
+  !                                                                 !
+  ! UNIT  FILE      PARAMETER(S)    DATA REPRESENTATION             !
+  !                                                                 !
+  ! 15    fort.15   U,V,ETA,T,PS,                                   !
+  !                 Q,SD,MSL,TCC,                                   !
+  !                 10U,10V,2T,2D,  regular lamda phi grid          !
+  !                 LSP,CP,SSHF,                                    !
+  !                 SSR,EWSS,NSSS                                   !
+  !                                                                 !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !
+
+  USE PHTOGR
+  USE GRTOPH
+  USE FTRAFO
+  USE RWGRIB2
+  USE GRIB_API
+
+  IMPLICIT NONE
+
+  REAL, ALLOCATABLE, DIMENSION (:,:)   :: Z
+  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: T, UV 
+  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: DIV, ETA
+  REAL, ALLOCATABLE, DIMENSION (:,:)   :: DPSDL, DPSDM
+  REAL, ALLOCATABLE, DIMENSION (:,:,:) :: PS,DPSDT
+  REAL, ALLOCATABLE, DIMENSION (:)     :: WSAVE,H
+  REAL, ALLOCATABLE, DIMENSION (:)     :: BREITE, GBREITE
+  REAL, ALLOCATABLE, DIMENSION (:)     :: AK, BK, PV
+  INTEGER                              :: NPV
+
+  ! Arrays for Gaussian grid calculations
+
+  REAL  :: X1,X2,RMS,MW,SIG,LAM
+  REAL,ALLOCATABLE :: CUA(:,:,:),CVA(:,:,:)
+
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: P,PP !,P2
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: XMN,HILFUV
+  REAL, ALLOCATABLE, DIMENSION (:)   :: LNPMN,LNPMN2,LNPMN3
+  REAL, ALLOCATABLE, DIMENSION (:)   :: WEIGHT
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: UGVG
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: DG, ETAG
+  REAL, ALLOCATABLE, DIMENSION (:,:) :: GWSAVE
+  REAL, ALLOCATABLE, DIMENSION (:)   :: PSG,HILF
+
+  ! end arrays for Gaussian grid calculations
+
+  INTEGER, ALLOCATABLE, DIMENSION (:) :: MLAT
+  INTEGER, ALLOCATABLE :: GIFAX(:,:)
+
+  REAL, PARAMETER   :: PI=ACOS(-1.D0)
+
+  REAL COSB,DAK,DBK,P00
+  REAL URLAR8,JMIN1,LLLAR8,MAXBMIN1,PIR8,DCOSB
+
+  INTEGER I,J,K,L,IERR,M,MK,NGI,NGJ
+  INTEGER LUNIT,LUNIT_OUT
+
+  INTEGER MAXL, MAXB, MLEVEL, LEVOUT,LEVMIN,LEVMAX
+  INTEGER MGAUSS,MSMOOTH, MNAUF,META
+  INTEGER MDPDETA,METAPAR
+  REAL RLO0, RLO1, RLA0, RLA1
+  CHARACTER*300 MLEVELIST
+
+  INTEGER MAUF, MANF,IFAX(10)
+
+  INTEGER igrib,iret
+
+  CHARACTER*80 FILENAME
+
+  NAMELIST /NAMGEN/ &
+       MAXL, MAXB,  &
+       MLEVEL,MLEVELIST,MNAUF,METAPAR, &
+       RLO0, RLO1, RLA0, RLA1, &
+       MGAUSS,MSMOOTH,META,&
+       MDPDETA
+
+  read (4,NAMGEN)
+
+  PRINT*, 'MAXL= ', MAXL, ' RLO0 = ', RLO0, ' RLO1 = ', RLO1
+  MAUF=INT(360.*(REAL(MAXL)-1.)/(RLO1-RLO0)+0.0001)
+  PRINT*, 'MAUF= ', MAUF
+
+  MANF=INT(REAL(MAUF)/360.*(360.+RLO0)+1.0001)
+  PRINT*, 'MANF= ', MANF
+
+  IF(MANF .gt. MAUF) MANF=MANF-MAUF
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !     GAUSS STUFF !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  ! computation of etadot on gaussian grid
+
+! a new grib message is loaded from an existing sample.
+! we want a regular gaussian grid
+! set the environment variable GRIB_SAMPLES_PATH
+  write (filename, "(A13,I0,A3,I0)") "EI_regular_gg",MNAUF+1,"_ml", MLEVEL
+  call grib_new_from_samples(igrib, FILENAME)
+
+  call grib_get(igrib,'numberOfPointsAlongAMeridian', NGJ)
+  ALLOCATE (MLAT(NGJ))
+
+  !     get as a integer
+  call grib_get(igrib,'pl', MLAT)
+  print*,'Number of points along a meridian  = ', NGJ
+  NGI=SUM(MLAT)
+
+  call grib_get(igrib,'numberOfVerticalCoordinateValues',mk)
+  print*, 'numberOfVerticalCoordinateValues = ', mk
+
+  IF(mk/2-1 .ne. MLEVEL) THEN 
+     WRITE(*,*) 'FATAL: Number of model levels',mk, &
+          ' does not agree with', MLEVEL,' in namelist'
+     STOP
+  ENDIF
+
+  call grib_get_size(igrib,'pv',npv)
+  allocate(pv(npv))
+  call grib_get(igrib,'pv',pv)
+  ALLOCATE(AK(NPV/2))
+  ALLOCATE(BK(NPV/2))
+  AK=pv(:NPV/2)
+  BK=pv(NPV/2+1:)
+  deallocate(pv)
+! END GAUSS INFO
+
+
+
+  ! Initialization of Legendre transform on LAT/LON grid.
+
+  ALLOCATE (BREITE(MAXB))
+  ALLOCATE (Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB))
+  !$OMP PARALLEL DO
+  DO  J=1,MAXB
+     BREITE(J)=SIN((RLA1-(J-1.D0)*(RLA1-RLA0)/(MAXB-1))* PI/180.D0)
+     CALL PLGNFA(MNAUF,BREITE(J),Z(0,J))
+  ENDDO
+  !$OMP END PARALLEL DO
+
+
+  ! Initialisation of fields for  FFT and Legendre transformation
+  !	to  Gaussian grid and back to phase space
+  ALLOCATE (GBREITE(NGJ),WEIGHT(NGJ))
+  ALLOCATE (P(0:((MNAUF+3)*(MNAUF+4))/2,NGJ/2))
+  ALLOCATE (PP(NGJ/2,0:((MNAUF+3)*(MNAUF+4))/2))
+  X1=-1.D0
+  X2=1.D0
+  CALL GAULEG(X1,X2,GBREITE,WEIGHT,NGJ)
+
+  !$OMP PARALLEL DO PRIVATE(M)
+  DO J=1,NGJ/2
+     CALL PLGNFA(MNAUF,GBREITE(J),P(:,J))
+     DO M=0,(MNAUF+3)*(MNAUF+4)/2
+        PP(J,M)=P(M,J)
+     ENDDO
+  ENDDO
+  !$OMP END PARALLEL DO
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!   READ LNSP
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  FILENAME='fort.12'
+  ALLOCATE (LNPMN(0:(MNAUF+1)*(MNAUF+2)-1))
+  CALL READSPECTRAL(FILENAME,LNPMN,MNAUF,1,MLEVEL,(/152/))
+  ! Call set99 an initialization routine that must be called once
+  ! before a sequence of calls to the fft routines
+  ALLOCATE (WSAVE(4*MAUF+15))
+  ALLOCATE (PS(MAXL, MAXB,1))
+  CALL SET99(WSAVE,IFAX,mauf)
+  CALL PHGCUT(LNPMN,PS,WSAVE,IFAX,Z, &
+       MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
+  CALL STATIS(MAXL,MAXB,1,EXP(PS),RMS,MW,SIG)
+  WRITE(*,'(A12,3F12.4)') 'STATISTICS : ',RMS,MW,SIG
+
+  ALLOCATE (GWSAVE(8*NGJ+15,NGJ/2))
+  ALLOCATE (GIFAX(10,NGJ))
+  DO J=1,NGJ/2
+     CALL SET99(GWSAVE(1,J),GIFAX(1,J),MLAT(J))
+  ENDDO
+  ALLOCATE (PSG(NGI),HILF(NGI))
+  ALLOCATE (LNPMN2(0:(MNAUF+1)*(MNAUF+2)-1))
+  CALL PHGR213(LNPMN,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
+  PSG=HILF
+  CALL GRPH213(LNPMN2,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT, &
+       MNAUF,NGI,NGJ,1)
+  CALL PHGR213(LNPMN2,HILF,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,1)
+
+
+  HILF=exp(PSG)-exp(HILF)
+
+  CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
+  WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
+
+  PSG=EXP(PSG)
+  HILF=PSG
+  CALL STATIS(NGI,1,1,HILF,RMS,MW,SIG)
+  WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!   READ U/V
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  FILENAME='fort.10'
+  ALLOCATE (XMN(0:(MNAUF+1)*(MNAUF+2)-1, 2*MLEVEL))
+  CALL READSPECTRAL(FILENAME, &
+       XMN,MNAUF,2*MLEVEL,MLEVEL,(/131,132/)) 
+  ! Transform winds on gaussian grid
+  ALLOCATE (UGVG(NGI, 2*MLEVEL))
+  CALL PHGR213(XMN,UGVG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,2*MLEVEL)
+  ALLOCATE (CUA(2,4,MLEVEL))
+  ALLOCATE (CVA(2,4,MLEVEL))
+  DO K=1,MLEVEL
+     ! North Pole
+     CALL JSPPOLE(XMN(:,K),1,MNAUF,.TRUE.,CUA(:,:,K))
+     CALL JSPPOLE(XMN(:,MLEVEL+K),1,MNAUF,.TRUE.,CVA(:,:,K))
+     ! South Pole
+     CALL JSPPOLE(XMN(:,K),-1,MNAUF,.TRUE.,CUA(:,3:4,K))
+     CALL JSPPOLE(XMN(:,MLEVEL+K),-1,MNAUF,.TRUE.,CVA(:,3:4,K))
+  ENDDO
+
+  IF(MSMOOTH .ne. 0) THEN
+    DO K=1,2*MLEVEL
+       CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
+    ENDDO
+  ENDIF
+  ALLOCATE (UV(MAXL, MAXB, 2*MLEVEL))
+  CALL PHGCUT(XMN,UV,WSAVE,IFAX,Z, &
+       MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,2*MLEVEL)
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!   READ Divergence
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  FILENAME='fort.13'
+  CALL READSPECTRAL(FILENAME,XMN,MNAUF,MLEVEL,MLEVEL,(/155/))
+  ! Tranform horizontal divergence on gaussian grid
+  ALLOCATE (DG(NGI,MLEVEL))
+  CALL PHGR213(XMN,DG,GWSAVE,GIFAX,P,MLAT,MNAUF,NGI,NGJ,MLEVEL)
+
+
+  ! compute gradient of LNSP (log of surface pressure) on gaussian grid 
+  ALLOCATE (H(0:(MNAUF+2)*(MNAUF+3)/2))
+  ALLOCATE (DPSDL(NGI,1),DPSDM(NGI,1))
+  CALL PHGRAD(LNPMN,DPSDL,DPSDM,GWSAVE,GIFAX,P,H,MLAT,MNAUF,NGI,NGJ,1)
+
+  ! Compute vertical weed spind on gaussian grid
+  ALLOCATE (ETAG(NGI,MLEVEL))
+  CALL CONTGL(HILF,DPSDL,DPSDM,DG,UGVG(:,1),UGVG(:,MLEVEL+1), &
+       GBREITE,ETAG,MLAT,AK,BK,NGI,NGJ,MLEVEL)
+
+
+  CALL GRPH213(XMN,ETAG,GWSAVE,GIFAX,PP,WEIGHT,MLAT, &
+       MNAUF,NGI,NGJ,MLEVEL)
+  IF(MSMOOTH .ne. 0) THEN
+    DO K=1,MLEVEL
+       CALL SPFILTER(XMN(:,K),MNAUF,MSMOOTH)
+    ENDDO
+  ENDIF
+  ALLOCATE (ETA(MAXL,MAXB,MLEVEL))
+  CALL PHGCUT(XMN,ETA,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
+
+  CALL GRPH213(XMN,HILF,GWSAVE,GIFAX,PP,WEIGHT,MLAT, MNAUF,NGI,NGJ,1)
+
+  IF(MSMOOTH .ne. 0) CALL SPFILTER(XMN(:,1),MNAUF,MSMOOTH)
+  ALLOCATE (DPSDT(MAXL, MAXB,1))
+  CALL PHGCUT(XMN,DPSDT,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
+
+  CALL STATIS(MAXL,MAXB,1,DPSDT,RMS,MW,SIG)
+  WRITE(*,'(A12,3F11.4)') 'STATISTICS DPSDT: ',RMS,MW,SIG
+
+  CALL GRPH213(XMN,PSG,GWSAVE,GIFAX,PP,WEIGHT,MLAT,MNAUF,NGI,NGJ,1)
+  CALL PHGCUT(XMN,PS,WSAVE,IFAX,Z,MNAUF,MNAUF,MAUF,MANF,MAXL,MAXB,1)
+
+  CALL STATIS(MAXL,MAXB,1,PS,RMS,MW,SIG)
+  WRITE(*,'(A12,3F11.4)') 'STATISTICS: ',RMS,MW,SIG
+
+114 DEALLOCATE(HILF,PSG,DPSDL,DPSDM,ETAG,DG,LNPMN)
+
+  DEALLOCATE(PP,P,UGVG,MLAT,GBREITE,WEIGHT,GWSAVE,XMN)
+
+  ! CREATE FILE VERTICAL.EC NEEDED BY POP MODEL
+
+  open(21,file='VERTICAL.EC')
+  write(21,'(a)')
+  write(21,'(a)') 'VERTICAL DISCRETIZATION OF POP MODEL'
+  write(21,'(a)')
+  write(21,'(i3,a)') MLEVEL,'   number of layers'
+  write(21,'(a)')
+  write(21,'(a)') '* A(NLEV+1)'
+  write(21,'(a)')
+  do  i=1,MLEVEL+1
+     write(21,'(f18.12)') AK(I)
+  enddo
+  write(21,'(a)')
+  write(21,'(a)') '* B(NLEV+1)'
+  write(21,'(a)')
+  do  i=1,MLEVEL+1
+     write(21,'(f18.12)') BK(I)
+  enddo
+  close(21)
+
+
+  !     Calculation of etadot in CONTGL needed scaled winds (ucosphi,vcosphi)
+  !     Now we are transforming back to the usual winds. 
+  ALLOCATE (HILFUV(2*MAXL,2))
+  DO K=1,MLEVEL
+     DO J=2,MAXB-1
+        COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
+        UV(:,J,K)=UV(:,J,K)/COSB
+        UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
+     ENDDO
+     ! special treatment for poles, if necessary. 
+     DO J=1,MAXB,MAXB-1
+        COSB=SQRT(1.0-(BREITE(J))*(BREITE(J)))
+        if(1.0-BREITE(J)*BREITE(J) .gt. 0 .OR. MGAUSS .NE. 1) then
+           IF(RLA0 .EQ. -90.0 .AND. J .EQ. MAXB .OR. &
+                RLA1 .EQ. 90.0 .AND. J .EQ. 1) then
+              UV(:,J,K)=UV(:,J,K)*1.D6
+              UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)*1.D6
+           else
+              UV(:,J,K)=UV(:,J,K)/COSB
+              UV(:,J,MLEVEL+K)=UV(:,J,MLEVEL+K)/COSB
+           endif
+        else
+           HILFUV(5:MAXL,:)=0.
+           HILFUV(1:2,:)=0.
+           IF(J.EQ.MAXB) THEN
+              ! Suedpol
+              HILFUV(3:4,1)=CUA(:,4,K)
+              HILFUV(3:4,2)=CVA(:,4,K)
+           ELSE
+              ! Nordpol
+              HILFUV(3:4,1)=CUA(:,2,K)
+              HILFUV(3:4,2)=CVA(:,2,K)
+           ENDIF
+           CALL RFOURTR(HILFUV(:,1),WSAVE,IFAX,MAXL/2-1,MAXL)
+           DO I=0,MAXL-1
+              IF(MANF+I.LE.MAXL) THEN
+                 UV(I+1,J,K)=HILFUV(MANF+I,1)
+              ELSE
+                 UV(I+1,J,K)=HILFUV(MANF-MAXL+I,1)
+              ENDIF
+           ENDDO
+           CALL RFOURTR(HILFUV(:,2),WSAVE,IFAX,MAXL/2-1,MAXL)
+           DO I=0,MAXL-1
+              IF(MANF+I.LE.MAXL) THEN
+                 UV(I+1,J,MLEVEL+K)=HILFUV(MANF+I,2)
+              ELSE
+                 UV(I+1,J,MLEVEL+K)=HILFUV(MANF-MAXL+I,2)
+              ENDIF
+           ENDDO
+        endif
+     ENDDO
+  ENDDO
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !                      READING OF T                      !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ALLOCATE (T(MAXL, MAXB, MLEVEL))
+  FILENAME='fort.11'
+  CALL READLATLON(FILENAME,T,MAXL,MAXB,MLEVEL,(/130/))
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !                    WRITE MODEL LEVEL DATA TO fort.15            !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! open output file
+  call grib_open_file(LUNIT,'fort.15','w')
+
+  ! we use temperature on lat/lon on model levels as template for model level data
+  LUNIT_OUT=0
+  call grib_open_file(LUNIT_OUT,'fort.11','r')
+  call grib_new_from_file(LUNIT_OUT,igrib, iret)
+  call grib_close_file(LUNIT_OUT)
+
+
+  CALL WRITELATLON(LUNIT,igrib,UV(:,:,1),MAXL,MAXB,MLEVEL,MLEVELIST,131)
+
+  CALL WRITELATLON(LUNIT,igrib,UV(:,:,MLEVEL+1),MAXL,MAXB,MLEVEL,MLEVELIST,132)
+
+  CALL WRITELATLON(LUNIT,igrib,ETA,MAXL,MAXB,MLEVEL,MLEVELIST,METAPAR)
+
+  CALL WRITELATLON(LUNIT,igrib,T,MAXL,MAXB,MLEVEL,MLEVELIST,130)
+
+  CALL WRITELATLON(LUNIT,igrib,PS,MAXL,MAXB,1,'1',134)
+
+  call grib_close_file(LUNIT)
+
+
+2000 STOP 'SUCCESSFULLY FINISHED CONVERT_PRE: CONGRATULATIONS'
+3000 STOP 'ROUTINE CONVERT_PRE: ERROR'
+9999 stop 'ROUTINE CONVERT_PRE: ERROR'
+
+CONTAINS
+
+  SUBROUTINE STATIS (NI,NJ,NK,PHI,RMS,MW,SIG)
+    IMPLICIT NONE
+
+    integer, intent(in) :: NI, NJ, NK
+
+    REAL PHI(NI,NJ,NK),SIG,MW,RMS,P
+    INTEGER N
+
+    N=NI*NJ*NK
+
+    RMS=0.
+    MW=0.
+
+    DO  I=1,NI
+       DO  J=1,NJ
+          DO  K=1,NK
+             P=PHI(I,J,K)
+             RMS=RMS+P*P
+             MW=MW+P
+          ENDDO
+       ENDDO
+    ENDDO
+
+    RMS=SQRT(RMS/N)
+    MW=MW/N
+
+    IF(RMS*RMS-MW*MW.LT.0.) THEN	  
+       SIG=0.0
+    ELSE
+       SIG=SQRT(RMS*RMS-MW*MW)
+    ENDIF
+
+    RETURN
+  END SUBROUTINE STATIS
+
+END PROGRAM PRECONVERT
