      PROGRAM OFFDIF
      USE MOD_ZA  ! HYCOM array I/O interface
      IMPLICIT NONE
C
C     WIND ARRAYS.
C
      INTEGER, ALLOCATABLE :: MSK(:,:)
      REAL*4,  ALLOCATABLE :: TXO(:,:),TYO(:,:),
     +                        TXN(:,:),TYN(:,:),
     +                        TXH(:,:),TYH(:,:)
C
      REAL*4  WDYNXT
C
C     NAMELISTS.
C
      CHARACTER*40     CTITLE
      REAL*4           SCALEN,SCALEO
      NAMELIST/WTITLE/ CTITLE
      NAMELIST/WOFSET/ SCALEN,SCALEO
C
C**********
C*
C  1) CREATE A WIND STRESS OFFSET, FROM A WIND SET'S ANNUAL CLIMATOLOGY
C      (OLD) AND THE DESIRED REPLACEMENT (NEW) ANNUAL CLIMATOLOGY.
C
C  2) REGION SPECIFICATION (FROM THE OCEAN MODEL):
C
C        IDM    = 1ST DIMENSION OF MAJOR MODEL ARRAYS
C        JDM    = 2ND DIMENSION OF MAJOR MODEL ARRAYS
C
C  3) INPUT:
C        ON UNIT  5:  NAMELISTS /WTITLE/ AND /WOFSET/
C        ON UNIT 44:  NEW .a/.b FORMAT MODEL MEAN TAUEWD
C        ON UNIT 45:  NEW .a/.b FORMAT MODEL MEAN TAUNWD
C        ON UNIT 54:  OLD .a/.b FORMAT MODEL MEAN TAUEWD
C        ON UNIT 55:  OLD .a/.b FORMAT MODEL MEAN TAUNWD
C     OUTPUT:
C        ON UNIT 64:  .a/.b FORMAT MODEL OFFSET TAUEWD
C        ON UNIT 65:  .a/.b FORMAT MODEL OFFSET TAUNWD
C
C     NAMELIST /WTITLE/:
C        CTITLE - 40 CHARACTER WIND OFFSET TITLE
C
C     NAMELIST /WOFSET/:
C        SCALEN - SCALE FACTOR FOR THE NEW ANNUAL CLIMATOLOGY, SEE (5)
C        SCALEO - SCALE FACTOR FOR THE OLD ANNUAL CLIMATOLOGY, SEE (5)
C
C  4) THE INPUT AND OUTPUT WIND STRESSES ARE ALWAYS ON THE HYCOM P-GRID
C      AND ARE IN HYCOM ARRAY FILE FORMAT.
C
C  5) 'SCALEN' AND 'SCALEO' ARE SUPPLIED TO SCALE THE DESIRED MEAN
C      SO THAT IT IS COMPATIBLE WITH THE WIND SET.  USUALLY FOR HYCOM,
C      SCALEN=SCALEO=1.0 (THE DEFAULT).  HOWEVER WHEN CONBINING
C      PARTIAL MEANS SCALEN-SCALEO=1.0 WOULD BE APPROPRIATE.
C      
C     THE OFFSET PRODUCED IS:  SCALEN*NEW - SCALEO*OLD.  
C
C     WHEN THIS IS ADDED TO THE WIND SET (ASSUMING SCALEO=1.0), 
C      IT WILL REPLACE THE OLD MEAN WITH SCALEN*NEW.
C
C  6) ALAN J. WALLCRAFT,  NAVAL RESEARCH LABORATORY,  JANUARY 2002.
C*
C**********
C
      INTEGER I,J
      REAL*4  XMIN,XMAX,XAVE,XRMS,YMIN,YMAX,YAVE,YRMS
      REAL*4  HMINA,HMINB,HMAXA,HMAXB
C
      CHARACTER PREAMBL(5)*79,CTITLO*40,CLINE*80
C
      CALL XCSPMD  !define idm,jdm
      ALLOCATE( TXO(IDM,JDM),TYO(IDM,JDM) )
      ALLOCATE( TXN(IDM,JDM),TYN(IDM,JDM) )
      ALLOCATE( TXH(IDM,JDM),TYH(IDM,JDM) )
C
C     READ IN, AND OUTPUT, THE WIND TITLE(S).
C
      CALL ZAIOST
C
      CALL ZHOPEN(6, 'FORMATTED', 'UNKNOWN', 0)
C
      CALL ZAIOPN('OLD', 44)
      CALL ZAIOPN('OLD', 45)
      CALL ZHOPEN(44, 'FORMATTED', 'OLD', 0)
      CALL ZHOPEN(45, 'FORMATTED', 'OLD', 0)
      READ(44,4101) PREAMBL
      READ(45,4101) PREAMBL
      CTITLO = PREAMBL(1)
      WRITE(6,6000) 'NEW MEAN:',CTITLO
      CALL ZHFLSH(6)
C
      CALL ZAIOPN('OLD', 54)
      CALL ZAIOPN('OLD', 55)
      CALL ZHOPEN(54, 'FORMATTED', 'OLD', 0)
      CALL ZHOPEN(55, 'FORMATTED', 'OLD', 0)
      READ(54,4101) PREAMBL
      READ(55,4101) PREAMBL
      CTITLO = PREAMBL(1)
      WRITE(6,6000) 'OLD MEAN:',CTITLO
      CALL ZHFLSH(6)
C
      READ( 5,WTITLE)
      WRITE(6,6000) 'OFFSET:  ',CTITLE
      WRITE(6,*)
      CALL ZHFLSH(6)
C
      CALL ZAIOPN('NEW', 64)
      CALL ZAIOPN('NEW', 65)
      CALL ZHOPEN(64, 'FORMATTED', 'NEW', 0)
      CALL ZHOPEN(65, 'FORMATTED', 'NEW', 0)
      PREAMBL(1) = CTITLE
      WRITE(64,4101) PREAMBL
      WRITE(65,4101) PREAMBL
      WRITE(6,*)
      WRITE(6, 4101) PREAMBL
      WRITE(6,*)
      CALL ZHFLSH(6)
C
C     NAMELIST /WOFSET/.
C
      SCALEO = 1.0
      SCALEN = 1.0
      READ( 5,WOFSET)
      WRITE(6,WOFSET)
      WRITE(6,*)
      CALL ZHFLSH(6)
C
C     READ IN THE NEW MEAN.
C
      READ(44,'(A)') CLINE
      I = INDEX(CLINE,'=')
      J = INDEX(CLINE,'month')
      IF     (J.NE.0) THEN
        READ(CLINE(I+1:),*) J,HMINB,HMAXB
      ELSE
        READ(CLINE(I+1:),*)   HMINB,HMAXB
      ENDIF
      CALL ZAIORD(TXN,MSK,.FALSE., HMINA,HMAXA, 44)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (txo):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
      CLOSE(UNIT=44)
      CALL ZAIOCL(44)
C
      READ(45,'(A)') CLINE
      I = INDEX(CLINE,'=')
      J = INDEX(CLINE,'month')
      IF     (J.NE.0) THEN
        READ(CLINE(I+1:),*) J,HMINB,HMAXB
      ELSE
        READ(CLINE(I+1:),*)   HMINB,HMAXB
      ENDIF
      CALL ZAIORD(TYN,MSK,.FALSE., HMINA,HMAXA, 45)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (tyo):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
      CLOSE(UNIT=45)
      CALL ZAIOCL(45)
C
      CALL MINMAX(TXN,IDM,JDM, XMIN,XMAX)
      CALL AVERMS(TXN,IDM,JDM, XAVE,XRMS)
      CALL MINMAX(TYN,IDM,JDM, YMIN,YMAX)
      CALL AVERMS(TYN,IDM,JDM, YAVE,YRMS)
      WRITE(6,*)
      WRITE(6,*) 'NEW MEAN STATISTICS.'
      WRITE(6,*)
      WRITE(6,8100) SCALEN*XMIN,SCALEN*XMAX,SCALEN*XAVE,SCALEN*XRMS
      WRITE(6,8200) SCALEN*YMIN,SCALEN*YMAX,SCALEN*YAVE,SCALEN*YRMS
      CALL ZHFLSH(6)
C
C     READ IN THE OLD MEAN.
C
      READ(54,'(A)') CLINE
      I = INDEX(CLINE,'=')
      J = INDEX(CLINE,'month')
      IF     (J.NE.0) THEN
        READ(CLINE(I+1:),*) J,HMINB,HMAXB
      ELSE
        READ(CLINE(I+1:),*)   HMINB,HMAXB
      ENDIF
      CALL ZAIORD(TXO,MSK,.FALSE., HMINA,HMAXA, 54)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (txo):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
      CLOSE(UNIT=54)
      CALL ZAIOCL(54)
C
      READ(55,'(A)') CLINE
      I = INDEX(CLINE,'=')
      J = INDEX(CLINE,'month')
      IF     (J.NE.0) THEN
        READ(CLINE(I+1:),*) J,HMINB,HMAXB
      ELSE
        READ(CLINE(I+1:),*)   HMINB,HMAXB
      ENDIF
      CALL ZAIORD(TYO,MSK,.FALSE., HMINA,HMAXA, 55)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (tyo):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
      CLOSE(UNIT=55)
      CALL ZAIOCL(55)
C
      CALL MINMAX(TXO,IDM,JDM, XMIN,XMAX)
      CALL AVERMS(TXO,IDM,JDM, XAVE,XRMS)
      CALL MINMAX(TYO,IDM,JDM, YMIN,YMAX)
      CALL AVERMS(TYO,IDM,JDM, YAVE,YRMS)
      WRITE(6,*)
      WRITE(6,*) 'OLD MEAN STATISTICS.'
      WRITE(6,*)
      WRITE(6,8100) XMIN,XMAX,XAVE,XRMS
      WRITE(6,8200) YMIN,YMAX,YAVE,YRMS
      WRITE(6,*)
      WRITE(6,*)
      CALL ZHFLSH(6)
C
C     OFFSET.
C
      DO 110 J= 1,JDM
        DO 111 I= 1,IDM
          TXH(I,J) = SCALEN*TXN(I,J) - SCALEO*TXO(I,J)
          TYH(I,J) = SCALEN*TYN(I,J) - SCALEO*TYO(I,J)
  111   CONTINUE
  110 CONTINUE
      CALL MINMAX(TXH,IDM,JDM, XMIN,XMAX)
      CALL AVERMS(TXH,IDM,JDM, XAVE,XRMS)
      CALL MINMAX(TYH,IDM,JDM, YMIN,YMAX)
      CALL AVERMS(TYH,IDM,JDM, YAVE,YRMS)
      WRITE(6,*)
      WRITE(6,*) 'OFFSET STATISTICS.'
      WRITE(6,*)
      WRITE(6,8100) XMIN,XMAX,XAVE,XRMS
      WRITE(6,8200) YMIN,YMAX,YAVE,YRMS
      WRITE(6,*)
      WRITE(6,*)
C
C     WRITE OUT OFFSET.
C
      CALL ZAIOWR(TXH,MSK,.FALSE., HMINA,HMAXA, 64, .FALSE.)
      WRITE(64,4102) ' tau_ewd',1,HMINA,HMAXA
      WRITE( 6,4102) ' tau_ewd',1,HMINA,HMAXA
      CALL ZAIOCL(64)
      CLOSE(UNIT=64)
C
      CALL ZAIOWR(TYH,MSK,.FALSE., HMINA,HMAXA, 65, .FALSE.)
      WRITE(65,4102) ' tau_nwd',1,HMINA,HMAXA
      WRITE( 6,4102) ' tau_nwd',1,HMINA,HMAXA
      CALL ZAIOCL(65)
      CLOSE(UNIT=65)
C
      WRITE(6,*)
      CALL ZHFLSH(6)
      STOP
C
 4101 FORMAT(A79)
 4102 FORMAT(A,': month,range = ',I2.2,1P2E16.7)
 6000 FORMAT(1X,A,2X,A40 //)
 6300 FORMAT(10X,'WRITING OFFSET RECORD' /)
 8100 FORMAT(1X,'TAU-X:  MIN =',F10.4,'   MAX =',F10.4,
     +              '    AVE =',F10.4,'   RMS =',F10.4)
 8200 FORMAT(1X,'TAU-Y:  MIN =',F10.4,'   MAX =',F10.4,
     +              '    AVE =',F10.4,'   RMS =',F10.4)
C     END OF PROGRAM OFFDIF.
      END
