C C Author: Ken w/ Sonali for the main procedure, but the HNTRP code C comes from Gary via email from Reto via email from Mark C C Purpose: Convert PRISM data to Model E grid resolution C C Method: C 1) Get PRISM Excel file C 2) Manually save each month as TXT file C 3) Convert to UNIX line endings: for v in *.txt; do mac2unix $v; done C 4) Now comes this code (read it, comments embedded) C 5) Use output in Model E! C C Run it once and move p2e.dat to the SST or SICE file. Run it a 2nd C time w or w/or SICE set and move p2e.dat to the appropriate file. C PROGRAM PRISM2E CHARACTER*3 MONTH(12)/'jan','feb','mar','apr','may', * 'jun','jul','aug','sep','oct','nov','dec'/ INTEGER IM,JM,I,J PARAMETER (IM=179, JM=90) REAL*4 PDATA(IM,JM), EDATA(72,46), WTA(IM,JM), PLINE(IM) INTEGER DLINE(IM+2), DIPOINT C "D" is "dummy" as in we will discard it. CHARACTER*80 HEADER INTEGER SICE COMMON SICE C Sea Ice? Change from 0 to 1 SICE = 1 WRITE(*,*) "Processing as SICE not SST" C Call the setup routine C$$$Input: IMA,JMA = array dimensions for A grid C$$$ OFFIA = number of grid boxes from IDL to left edge of gridbox I=1 C$$$ DIVJA = number of whole grid boxes from SP to NP C$$$ IMB,JMB = array dimensions for B grid C$$$ OFFIB = number of grid boxes from IDL to left edge of gridbox I=1 C$$$ DIVJB = number of whole grid boxes from SP to NP C$$$ SKIP = value to be put in B array when a B grid box is a subset C$$$ of A grid boxes with WTA = 0 IMA = 179 JMA = 90 OFFIA = 0 DIVJA = 90 IMB = 72 JMB = 46 OFFIB = 0 DIVJB = 46 SKIP = -999 CALL HNTRP0 (IMA,JMA,OFFIA,DIVJA,IMB,JMB,OFFIB,DIVJB,SKIP) HEADER='PRISM2 SST YEAR 4x5 conv by Mankoff. Mean=End-of-Month.' WRITE(*,*) "Converting PRISM to E" DO 10 M=1,12 C Open the data file C WRITE(*,*) "Reading: "//MONTH(M)//".txt" OPEN (100, FILE=TRIM(MONTH(M))//".txt", STATUS='OLD') C read the 1st line and discard READ(100,FMT=*) DLINE C read each line. Throw out 1st value (lat). Store the rest in 2D array DO 11 J=1,90 READ(100,FMT=*) DIPOINT, PLINE DO 12 I=1,179 IF (SICE.EQ.1) THEN IF (PLINE(I).EQ.-1.8) THEN PLINE(I)=1 ELSE PLINE(I)=0 ENDIF ENDIF PDATA(I,J) = PLINE(I) WTA(I,J) = 1 IF (PLINE(I).EQ.-999.) WTA(I,J) = 0 12 CONTINUE 11 CONTINUE CLOSE (100) C WRITE(*,*) "Converting : "//MONTH(M)//" PRISM to E" C Do the conversion CALL HNTRPP(WTA,PDATA,EDATA) C write out file IF (M.EQ.1) OPEN(UNIT=101, FILE='./p2e.dat', * FORM='unformatted') WRITE(101) HEADER, EDATA, EDATA C write out edata twice, because the file contains month and eom arrays 10 ENDDO CLOSE(101) CALL EOM END SUBROUTINE EOM REAL*4 EDATA12(72,46,12), EDATA0(72,46), EDATA012(72,46,12) REAL*4 EDATA(72,46), DPICEX(72,46) CHARACTER*80 HEADER, HEADERS(12), HEADERI INTEGER SICE COMMON SICE WRITE(*,*) "Converting 2nd array to EOM data" HEADERI='Factor needed to relate ice depth to ice cover' C addition: Re-open file, read in everything, and C calculate End Of Month. OPEN(UNIT=102, FILE='./p2e.dat',FORM='unformatted') DO M=1,12 READ(102) HEADER, EDATA, EDATA0 DO I=1,IM*JM HEADERS(M) = HEADER EDATA12(I,1,M) = EDATA(I,1) EDATA012(I,1,M) = EDATA0(I,1) END DO END DO CLOSE(102) C EDATA12 and EDATA012 now contain 12 months of data C Taken from OST4X5.B.1958-65avg.f MM1 = 10 M = 11 MP1 = 12 DO 90 MP2=1,12 DO 80 I=1,IM*JM EDATA012(i,1,m)=(7.*(EDATA12(i,1,m)+EDATA12(i,1,mp1)) * -(EDATA12(i,1,mm1)+EDATA12(i,1,mp2)))/12. IF(EDATA012(i,1,m).lt.-1.8) EDATA012(i,1,m)=-1.8 IF(EDATA12(i,1,m).le.-1.79) EDATA012(i,1,m)=-1.8 IF(EDATA12(i,1,mp1).lt.-1.79) EDATA012(i,1,m)=-1.8 80 CONTINUE MM1 = M M = MP1 90 MP1 = MP2 C End of OST4X5 code C Re-write SST file OPEN(UNIT=103, FILE='./p2e.dat',FORM='unformatted') C find ice depth. From SICE4X5.B.1951-65avg.f do i=1,im*jm DPICEX(i,1)=0. end do do m=1,12 do i=1,im*jm if(EDATA12(i,1,m).gt.0.) DPICEX(i,1)=DPICEX(i,1)+1. end do end do do i=1,im*jm if(DPICEX(i,1).gt.11.5) DPICEX(i,1)=16 DPICEX(i,1)=.25*sqrt(DPICEX(i,1)) end do C ICE DEPTH CALCULATIONS DONE C calculate End Of Month. DO M=1,12 DO I=1,IM*JM EDATA(I,1) = EDATA12(I,1,M) EDATA0(I,1) = EDATA012(I,1,M) END DO HEADER = HEADERS(M) IF ((SICE.EQ.1).AND.(M.EQ.1)) THEN WRITE(*,*) "Writing Ice Depth First" WRITE(103) HEADERI,DPICEX END IF WRITE(103) HEADER, EDATA, EDATA0 END DO CLOSE(103) END C$$$ SUBROUTINE P2E C$$$ implicit none C$$$ INTEGER*4 IMA,JMA,OFFIA,DIVJA,IMB,JMB,OFFIB,DIVJB,SKIP C$$$ INTEGER*4 T C$$$ REAL*4 PDATA(179,90), EDATA(72,46) C$$$ CHARACTER*80 HEADER C$$$ IMA = 179 C$$$ JMA = 90 C$$$ OFFIA = 0 C$$$ DIVJA = 90 C$$$ IMB = 72 C$$$ JMB = 46 C$$$ OFFIB = 0 C$$$ DIVJB = 46 C$$$ SKIP = -999 C$$$ CALL HNTRP0 (IMA,JMA,OFFIA,DIVJA, IMB,JMB,OFFIB,DIVJB, SKIP) C$$$ HEADER='PRISM2 SST YEAR 4x5 conv by Mankoff. Mean=End-of-Month.' C$$$C read in the file C$$$ OPEN(UNIT=100, FILE='pb.dat', STATUS='OLD',FORM='unformatted') C$$$ DO 42 T=1,12 C$$$ READ(100) PDATA C$$$C SISE is treated different than SST. Don't forget to change the file name! C$$$C DO 10 I=1,179 C$$$C DO 11 J=1,90 C$$$C IF (PDATA(I,J).EQ.-1.8) PDATA(I,J)=1.ELSE.PDATA(I,J)=0 C$$$C 11 CONTINUE C$$$C 10 CONTINUE C$$$C Do the conversion C$$$ CALL HNTRPP(1,PDATA,EDATA) C$$$C DO 11 J=1,1000 C$$$C WRITE(6,*) EDATA(0,J) C$$$C EDATA(0,J) = EDATA(0,J) C$$$C WRITE(6,*) J C$$$C 11 CONTINUE C$$$C write out file C$$$ IF (T.EQ.1) OPEN(UNIT=101, FILE='./p2eSST.dat',FORM='unformatted') C$$$ WRITE(101) HEADER, EDATA, EDATA C$$$C write out edata twice, because the file contains month and eom arrays C$$$ 42 CONTINUE C$$$ CLOSE(101) C$$$ END C**** HNTRPS.FOR Horizontal Interpolation Program REAL*8 2004/03/17 C**** SUBROUTINE HNTRP0 (INA,JNA,OFFIA,DIVJA, * INB,JNB,OFFIB,DIVJB, SKIB) C**** C**** HNTRP performs a horizontal interpolation of per unit area or per C**** unit mass quantities defined on grid A, calculating the quantity C**** on grid B. B grid values that cannot be calculated because the C**** covering A grid boxes have WTA = 0, are set to the value of SKIP. C**** The area weighted integral of the quantity is conserved. C**** IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (TWOPI=6.283185307179586477d0) C REAL*4 WTA(INA,JNA),A(INA,JNA),B(INB,JNB), REAL*4 WTA(*), A(*), B(*), * OFFIA,DIVJA, OFFIB,DIVJB, SKIB,SKIP COMMON /HNTRCB/ SINA(0:361),SINB(0:361), * FMIN(720),FMAX(720),GMIN(361),GMAX(361), * IMIN(720),IMAX(720),JMIN(362),JMAX(362) LOGICAL*4 QMPOLE DATA IMA,JMA,IMB,JMB/4*0/, SKIP/0/ C**** IMA = INA JMA = JNA IMB = INB JMB = JNB SKIP = SKIB IF(IMA.lt.1 .or. IMA.gt.720 .or. JMA.lt.1 .or. JMA.gt.361 .or. * IMB.lt.1 .or. IMB.gt.720 .or. JMB.lt.1 .or. JMB.gt.361) * GO TO 400 C**** C**** Partitions in the I direction C**** RIA = longitude in degrees of right edge of grid box IA on grid A C**** RIB = longitude in degrees of right edge of grid box IB of grid B C**** IMIN(IB) = box on grid A containing left edge of box IB on B C**** IMAX(IB) = box on grid A containing right edge of box IB on B C**** FMIN(IB) = fraction of box IMIN(IB) on A that is left of box IB C**** FMAX(IB) = fraction of box IMAX(IB) on A that is right of box IB C**** DIA = 360d0/IMA DIB = 360d0/IMB IA = 1 RIA = (IA+OFFIA)*DIA - 360 IB = IMB DO 150 IBP1=1,IMB RIB = (IBP1-1+OFFIB)*DIB 110 IF(RIA-RIB) 120,130,140 120 IA = IA + 1 RIA = RIA + DIA GO TO 110 C**** Right edge of A box IA and right edge of B box IB coincide 130 IMAX(IB) = IA FMAX(IB) = 0. IA = IA + 1 RIA = RIA + DIA IMIN(IBP1) = IA FMIN(IBP1) = 0. GO TO 150 C**** A box IA contains right edge of B box IB 140 IMAX(IB) = IA FMAX(IB) = (RIA-RIB)/DIA IMIN(IBP1) = IA FMIN(IBP1) = 1.-FMAX(IB) 150 IB = IBP1 IMAX(IMB) = IMAX(IMB) + IMA C WRITE (0,*) 'IMIN=',(IMIN(I),I=1,IMB) C WRITE (0,*) 'IMAX=',(IMAX(I),I=1,IMB) C WRITE (0,*) 'FMIN=',(FMIN(I),I=1,IMB) C WRITE (0,*) 'FMAX=',(FMAX(I),I=1,IMB) C**** C**** Partitions in the J direction C**** C**** RJA = latitude in radians at top edge of box JA on grid A C**** SINA(JA) = sine of latitude of top edge of box JA on grid A OFFJA = (DIVJA-JMA)/2. DJA = .5*TWOPI/DIVJA DO 210 JA=1,JMA-1 RJA = (JA+OFFJA)*DJA - .25*TWOPI 210 SINA(JA) = DSIN(RJA) SINA(0) = -1. SINA(JMA)= 1. C**** RJB = latitude in radians at top edge of box JB on grid B C**** SINB(JB) = sine of latitude of top edge of box JB on grid B OFFJB = (DIVJB-JMB)/2. DJB = .5*TWOPI/DIVJB DO 220 JB=1,JMB-1 RJB = (JB+OFFJB)*DJB - .25*TWOPI 220 SINB(JB) = DSIN(RJB) SINB(0) = -1. SINB(JMB)= 1. C**** C**** JMIN(JB) = index of box of A that contains bottom edge of box JB C**** JMAX(JB) = index of box of A that contains top edge of box JB C**** GMIN(JB) = fraction of box JMIN(JB) on A grid that is below box JB C**** GMAX(JB) = fraction of box JMAX(JB) on A grid that is above box JB C**** JMIN(1) = 1 GMIN(1) = 0. JA = 1 DO 350 JB=1,JMB-1 310 IF(SINA(JA)-SINB(JB)) 320,330,340 320 JA = JA + 1 GO TO 310 C**** Top edge of A box JA and top edge of B box JB coincide 330 JMAX(JB) = JA GMAX(JB) = 0. JA = JA + 1 JMIN(JB+1) = JA GMIN(JB+1) = 0. GO TO 350 C**** A box JA contains top edge of B box JB 340 JMAX(JB) = JA GMAX(JB) = SINA(JA)-SINB(JB) JMIN(JB+1) = JA GMIN(JB+1) = SINB(JB)-SINA(JA-1) 350 CONTINUE JMAX(JMB) = JMA GMAX(JMB) = 0. C WRITE (0,*) 'JMIN=',(JMIN(J),J=1,JMB) C WRITE (0,*) 'JMAX=',(JMAX(J),J=1,JMB) C WRITE (0,*) 'GMIN=',(GMIN(J),J=1,JMB) C WRITE (0,*) 'GMAX=',(GMAX(J),J=1,JMB) RETURN C**** C**** Invalid parameters or dimensions out of range C**** 400 WRITE (0,940) IMA,JMA,OFFIA,DIVJA, IMB,JMB,OFFIB,DIVJB, SKIP STOP 400 940 FORMAT ('0Arguments received by HNTRP0 in order:'/ * 2I12,' = IMA,JMA = array dimensions for A grid'/ * E24.8,' = OFFIA = fractional number of grid boxes from', * ' IDL to left edge of grid box I=1'/ * E24.8,' = DIVJA = number of whole grid boxes from SP to NP'/ * 2I12,' = IMB,JMB = array dimensions for B grid'/ * E24.8,' = OFFIB = fractional number of grid boxes from', * ' IDL to left edge of grid box I=1'/ * E24.8,' = DIVJB = number of whole grid boxes from SP to NP'/ * E24.8,' = SKIP = value to be put in B array when B', * ' grid box is subset of A grid boxes with WTA = 0'/ * '0These arguments are invalid or out of range.') C**** ENTRY HNTRP (WTA,A,B) C**** C**** HNTRP performs the horizontal interpolation C**** Input: WTA = weighting array for values on the A grid C**** A = per unit area or per unit mass quantity C**** Output: B = horizontally interpolated quantity on B grid C**** QMPOLE = .FALSE. GO TO 500 ENTRY HNTRPP (WTA,A,B) C**** C**** HNTRPP is similar to HNTRP but polar values are replaced by C**** their longitudinal mean C**** QMPOLE = .TRUE. C**** C**** Interpolate the A grid onto the B grid C**** 500 DO 520 JB=1,JMB JAMIN = JMIN(JB) JAMAX = JMAX(JB) DO 520 IB=1,IMB IJB = IB + IMB*(JB-1) WEIGHT= 0. VALUE = 0. IAMIN = IMIN(IB) IAMAX = IMAX(IB) DO 510 JA=JAMIN,JAMAX G = SINA(JA)-SINA(JA-1) IF(JA.eq.JAMIN) G = G - GMIN(JB) IF(JA.eq.JAMAX) G = G - GMAX(JB) DO 510 IAREV=IAMIN,IAMAX IA = 1+MOD(IAREV-1,IMA) IJA = IA + IMA*(JA-1) F = 1. IF(IAREV.eq.IAMIN) F = F - FMIN(IB) IF(IAREV.eq.IAMAX) F = F - FMAX(IB) WEIGHT = WEIGHT + F*G*WTA(IJA) 510 VALUE = VALUE + F*G*WTA(IJA)*A(IJA) B(IJB) = SKIP IF(WEIGHT.ne.0.) B(IJB) = VALUE/WEIGHT 520 continue C**** C**** Replace individual values near the poles by longitudinal mean C**** IF(.NOT.QMPOLE) RETURN DO 630 JB=1,JMB,JMB-1 WEIGHT = 0. VALUE = 0. DO 610 IB=1,IMB IJB = IB + IMB*(JB-1) IF(B(IJB).eq.SKIP) GO TO 610 WEIGHT = WEIGHT + 1. VALUE = VALUE + B(IJB) 610 continue BMEAN = SKIP IF(WEIGHT.ne.0.) BMEAN = VALUE/WEIGHT DO 620 IB=1,IMB IJB = IB + IMB*(JB-1) 620 B(IJB) = BMEAN 630 continue RETURN END