subroutine INITLZ & (im,jm,platn,plonn,xt,yt,xm,ym,dxeta0,dyeta0,GLAT) !================================================================ ! *** This routine is the driver of all initialization packages !--------------------------------------------------------------------- REAL platn,plonn INTEGER IM,JM REAL,DIMENSION(IM,JM) :: DXU,DXV,DXT,DXM REAL,DIMENSION(IM,JM) :: DYU,DYV,DYT,DYM REAL,DIMENSION(IM,JM) :: GLAT,GLON REAL,DIMENSION(IM) :: xm,xt REAL,DIMENSION(JM) :: ym,yt REAL,DIMENSION(:),ALLOCATABLE :: XMN,XTN REAL,DIMENSION(:),ALLOCATABLE :: YMN,YTN REAL,DIMENSION(:,:),ALLOCATABLE :: FMAPU,FMAPV,FMAPT,FMAPM !---------------------------------------------------------------------- ! Initial startup !---------------------------------------------------------------------- PRINT*,' Initial start' ! Calculate grids, define topography and transform, and define ! initial fields. allocate(xmn(IM)) allocate(ymn(JM)) allocate(xtn(IM)) allocate(ytn(JM)) CALL GRIDSET & (1,IM,JM,platn,plonn,xmn,ymn,xtn,ytn) ! Fill latitude-longitude, map factor, and Coriolis arrays. CALL NEWGRID & (1,IM,JM,xmn,ymn,xtn,ytn,xt,yt,xm,ym) allocate(FMAPU(IM,JM)) allocate(FMAPV(IM,JM)) allocate(FMAPT(IM,JM)) allocate(FMAPM(IM,JM)) CALL POLARST & (IM,JM,platn,plonn,xm,ym,xt,yt,GLAT,GLON,FMAPU,FMAPV,FMAPT,FMAPM) CALL GRDSPC & (IM,JM,FMAPU,FMAPV,FMAPT,FMAPM,XMN,YMN,XTN,YTN & ,DXU,DXV,DXT,DXM,DYU,DYV,DYT,DYM) ! define equi-distant grid dxeta0=0.0 do i=1,IM-1 dxeta0=dxeta0+1./DXM(i,1) end do dxeta0=dxeta0/(im-1) dyeta0=0.0 do j=1,JM-1 dyeta0=dyeta0+1./DYM(1,j) end do dyeta0=dyeta0/(jm-1) deallocate(FMAPU) deallocate(FMAPV) deallocate(FMAPT) deallocate(FMAPM) deallocate(xmn) deallocate(ymn) deallocate(xtn) deallocate(ytn) return end subroutine INITLZ subroutine GRIDSET & (NGRA,IM,JM,platn,plonn,xmn,ymn,xtn,ytn) REAL XMN(IM),XTN(IM) REAL YMN(JM),YTN(JM) REAL deltaxn,deltayn,platn,plonn INTEGER NGRA,IM,JM INTEGER JDIM CHARACTER(30) LLNAME !!!! DIMENSION ZMNVC(-1:NZPMAX+1) NAMELIST /RAMSIN_LATLON/ & DELTAX,DELTAY,POLELAT,POLELON,CENTLAT,CENTLON !---------------------------------------------------------- write(LLNAME,101) 101 format('ramsin_latlon.name') irams=11 CLOSE(irams) OPEN(UNIT=irams,FILE=LLNAME,FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(3,*)' irams OPEN UNIT ERROR IER=',IER READ(irams,RAMSIN_LATLON) CLOSE(irams,STATUS='KEEP') IF(JM.EQ.1) then JDIM=0 ELSE JDIM=1 END IF ! Fill platn and plonn arrays, the value of polelat and polelon for each ! grid. platn = polelat plonn = polelon ! FILL ALL DELTAXN AND DELTAYN VALUES, AND FIND NINEST AND NJNEST VALUES ! IF SET TO ZERO IN NAMELIST deltaxn = deltax deltayn = deltay call ll_xy(centlat,centlon,platn,plonn,centx1,centy1) xmn(1) = centx1 - 0.5 * float(IM-2) * deltaxn ymn(1) = centy1 - 0.5 * float(JM-2) * deltayn if (JM .eq. 1) ymn(1) = centy1 ! Calculate xmn and ymn for the coarse grid(s) for an initial start. ! Calculate xtn and ytn for the coarse grid(s) for an initial start. DO I = 2,IM XMN(I) = XMN(I-1) + DELTAXN XTN(I) = .5 * (XMN(I) + XMN(I-1)) ENDDO XTN(1) = 1.5 * XMN(1) - .5 * XMN(2) DO J = 2,JM YMN(J) = YMN(J-1) + DELTAYN YTN(J) = .5 * (YMN(J) + YMN(J-1)) ENDDO YTN(1) = 1.5 * YMN(1) - .5 * YMN(2) !======== use this in prep_corr for z-correlations ============== ! Calculate zmn for the coarse grid(s). ! NESTZ1=0 ! NESTZ2=0 ! DZRAT=1.03 ! DZMAX=1000. ! DELTAZ=0. ! IF (DELTAZ .EQ. 0.) THEN ! DO K = 1,NNZP ! ZMN(K) = ZZ(K) ! ENDDO ! ELSE ! ZMN(1) = 0. ! ZMN(2) = DELTAZ ! DO K = 3,NNZP ! DZR = DZRAT ! ZMN(K) = ZMN(K-1) ! + + MIN(DZR * (ZMN(K-1) - ZMN(K-2)),DZMAX) ! ENDDO ! DELTAZN = DELTAZ ! ENDIF ! Fill scratch array znmvc with ZMN for coarse grid(s) ! ZTOP = ZMN(NNZP-1) ! DO K = 1,NNZP ! ZMNVC(K) = ZMN(K) ! ENDDO ! ZMNVC(0) = -(ZMNVC(2) - ZMNVC(1)) ** 2 ! + / (ZMNVC(3) - ZMNVC(2)) ! ZMNVC(-1) = ZMNVC(0) - (ZMNVC(1) - ZMNVC(0)) ** 2 ! + / (ZMNVC(2) - ZMNVC(1)) ! ZMNVC(NNZP+1) = ZMNVC(NNZP) ! + + (ZMNVC(NNZP) - ZMNVC(NNZP-1)) ** 2 ! + / (ZMNVC(NNZP-1) - ZMNVC(NNZP-2)) ! ! Compute ZTN values for all grids by geometric interpolation. ! ! DO K = 1,NNZP ! DZRFM = SQRT(SQRT((ZMNVC(K+1) - ZMNVC(K)) ! + / (ZMNVC(K-1) - ZMNVC(K-2)))) ! ZTN(K) = ZMNVC(K-1) ! + + (ZMNVC(K) - ZMNVC(K-1)) / (1. + DZRFM) ! ENDDO ! ! Compute other arrays based on the vertical grid. ! ! DO K = 1,NNZP-1 ! DZMN(K) = 1. / (ZTN(K+1) - ZTN(K)) ! ENDDO ! DO K = 2,NNZP ! DZTN(K)=1./(ZMN(K) - ZMN(K-1)) ! ENDDO ! DO K = 2,NNZP-1 ! DZM2N(K) = 1. / (ZMN(K+1) - ZMN(K-1)) ! DZT2N(K) = 1. / (ZTN(K+1) - ZTN(K-1)) ! ENDDO ! DZMN(NNZP)=DZMN(NNZP-1) ! + *DZMN(NNZP-1)/DZMN(NNZP-2) ! DZTN(1)=DZTN(2)*DZTN(2)/DZTN(3) ! DZM2N(1)=DZM2N(2) ! DZM2N(NNZP)=DZM2N(NNZP-1) ! DZT2N(1)=DZT2N(2) ! DZT2N(NNZP)=DZT2N(NNZP-1) ! ENDDO !================================================ return end subroutine GRIDSET subroutine POLARST & (N2,N3,platn,plonn,xm,ym,xt,yt,GLAT,GLON,FMAPU,FMAPV,FMAPT,FMAPM) REAL erad,c1 PARAMETER (erad = 6367000.) REAL,DIMENSION(N2,N3) :: GLAT,GLON REAL,DIMENSION(N2,N3) :: FMAPU ,FMAPV ,FMAPT ,FMAPM REAL,DIMENSION(N2,N3) :: FMAPUI,FMAPVI,FMAPTI,FMAPMI REAL,DIMENSION(N2) :: XM,XT REAL,DIMENSION(N3) :: YM,YT REAL platn,plonn,xm2,xt2,ym2,yt2 INTEGER N2,N3 ! Calculates map factors and inverse map factors at u,v,t,m-points and ! geographical lat/lon at t-points for a given polar stereographic grid c1 = (2. * erad) ** 2 DO J = 1,N3 DO I = 1,N2 xm2 = xm(i) * xm(i) xt2 = xt(i) * xt(i) ym2 = ym(j) * ym(j) yt2 = yt(j) * yt(j) FMAPT(I,J) = 1. + (xt2 + yt2) / c1 FMAPU(I,J) = 1. + (xm2 + yt2) / c1 FMAPV(I,J) = 1. + (xt2 + ym2) / c1 FMAPM(I,J) = 1. + (xm2 + ym2) / c1 FMAPUI(I,J) = 1.0 / FMAPU(I,J) FMAPVI(I,J) = 1.0 / FMAPV(I,J) FMAPTI(I,J) = 1.0 / FMAPT(I,J) FMAPMI(I,J) = 1.0 / FMAPM(I,J) call xy_ll(GLAT(I,J),GLON(I,J),platn,plonn,XT(I),YT(J)) ENDDO ENDDO return end subroutine POLARST subroutine GRDSPC & (N2,N3,FMAPU,FMAPV,FMAPT,FMAPM,XMN,YMN,XTN,YTN & ,DXU,DXV,DXT,DXM,DYU,DYV,DYT,DYM) REAL,DIMENSION(N2,N3) :: FMAPU,FMAPV,FMAPT,FMAPM REAL,DIMENSION(N2,N3) :: DXU,DXV,DXT,DXM REAL,DIMENSION(N2,N3) :: DYU,DYV,DYT,DYM REAL,DIMENSION(N2) :: XMN,XTN REAL,DIMENSION(N3) :: YMN,YTN INTEGER N2,N3 !--------------------------------------------------- DO J=1,N3 DO I=1,N2-1 DXU(I,J)=FMAPU(I,J)/(XTN(I+1)-XTN(I)) DXM(I,J)=FMAPM(I,J)/(XTN(I+1)-XTN(I)) ENDDO DXU(N2,J)=DXU(N2-1,J)*FMAPU(N2,J)/FMAPU(N2-1,J) DXM(N2,J)=DXM(N2-1,J)*FMAPM(N2,J)/FMAPM(N2-1,J) DO I=2,N2 DXV(I,J)=FMAPV(I,J)/(XMN(I)-XMN(I-1)) DXT(I,J)=FMAPT(I,J)/(XMN(I)-XMN(I-1)) ENDDO DXV(1,J)=DXV(2,J)*FMAPV(1,J)/FMAPV(2,J) DXT(1,J)=DXT(2,J)*FMAPT(1,J)/FMAPT(2,J) ENDDO DO I=1,N2 DO J=1,N3-1 DYV(I,J)=FMAPV(I,J)/(YTN(J+1)-YTN(J)) DYM(I,J)=FMAPM(I,J)/(YTN(J+1)-YTN(J)) ENDDO DYV(I,N3)=DYV(I,N3-1)*FMAPV(I,N3)/FMAPV(I,N3-1) DYM(I,N3)=DYM(I,N3-1)*FMAPM(I,N3)/FMAPM(I,N3-1) DO J=2,N3 DYU(I,J)=FMAPU(I,J)/(YMN(J)-YMN(J-1)) DYT(I,J)=FMAPT(I,J)/(YMN(J)-YMN(J-1)) ENDDO DYU(I,1)=DYU(I,2)*FMAPU(I,1)/FMAPU(I,2) DYT(I,1)=DYT(I,2)*FMAPT(I,1)/FMAPT(I,2) ENDDO return end subroutine GRDSPC subroutine NEWGRID & (NGR,IM,JM,xmn,ymn,xtn,ytn,xt,yt,xm,ym) ! +---------------------------------------------------------------- ! ! Fill the single and 1D variables that the rest of the model ! ! uses from the nest arrays and change grid level in the I/O. ! +---------------------------------------------------------------- REAL XMN(IM),XTN(IM),XT(IM),XM(IM) REAL YMN(JM),YTN(JM),YT(JM),YM(JM) INTEGER NGR,IM,JM ! GRID SPACINGS DO I=1,IM XT(I)=XTN(I) XM(I)=XMN(I) ENDDO DO J=1,JM YT(J)=YTN(J) YM(J)=YMN(J) ENDDO !=============================== ! DELTAZ=DELTAZN ! DO K=1,NZP ! ZT(K)=ZTN(K) ! ZM(K)=ZMN(K) ! DZM(K)=DZMN(K) ! DZT(K)=DZTN(K) ! DZM2(K)=DZM2N(K) ! DZT2(K)=DZT2N(K) ! ENDDO !================================ return end subroutine NEWGRID subroutine ll_xy(qlat,qlon,polelat,polelon,x,y) parameter (erad=6.367e6,erad2=1.2734e7,pi180=3.14159265/180.) ! Evaluate sine and cosine of latitude and longitude of pole point p and ! input point q. sinplat = sin(polelat * pi180) cosplat = cos(polelat * pi180) sinplon = sin(polelon * pi180) cosplon = cos(polelon * pi180) sinqlat = sin(qlat * pi180) cosqlat = cos(qlat * pi180) sinqlon = sin(qlon * pi180) cosqlon = cos(qlon * pi180) ! Compute (x3,y3,z3) coordinates where the origin is the center of the earth, ! the z axis is the north pole, the x axis is the equator and prime ! meridian, and the y axis is the equator and 90 E. ! For the pole point, these are: x3p = erad * cosplat * cosplon y3p = erad * cosplat * sinplon z3p = erad * sinplat ! For the given lat,lon point, these are: z3q = erad * sinqlat x3q = erad * cosqlat * cosqlon y3q = erad * cosqlat * sinqlon ! Transform q point from (x3,y3,z3) coordinates in the above system to ! polar stereographic coordinates (x,y,z): xq=-sinplon*(x3q-x3p)+cosplon*(y3q-y3p) yq=cosplat*(z3q-z3p)-sinplat*( cosplon*(x3q-x3p)+sinplon*(y3q-y3p) ) zq=sinplat*(z3q-z3p)+cosplat*( cosplon*(x3q-x3p)+sinplon*(y3q-y3p) ) ! Parametric equation for line from antipodal point at (0,0,-2 erad) to ! point q has the following parameter (t) value on the polar stereographic ! plane: t = erad2 / (erad2 + zq) ! This gives the following x and y coordinates for the projection of point q ! onto the polar stereographic plane: x = xq * t y = yq * t return end subroutine ll_xy subroutine xy_ll(qlat,qlon,polelat,polelon,x,y) parameter (erad=6.367e6,erad2=1.2734e7,pi180=3.14159265/180.) ! Evaluate sine and cosine of latitude and longitude of pole point p. sinplat = sin(polelat * pi180) cosplat = cos(polelat * pi180) sinplon = sin(polelon * pi180) cosplon = cos(polelon * pi180) ! Compute (x3,y3,z3) coordinates of the pole point where the origin is the ! center of the earth, the z axis is the north pole, the x axis is the ! equator and prime meridian, and the y axis is the equator and 90 E. x3p = erad * cosplat * cosplon y3p = erad * cosplat * sinplon z3p = erad * sinplat ! Compute distance d from given point R on the polar stereographic plane ! to the pole point P: d = sqrt (x ** 2 + y ** 2) ! Compute angle QCP where C is the center of the Earth. This is twice ! angle QAP where A is the antipodal point. Angle QAP is the same as ! angle RAP: alpha = 2. * atan2(d,erad2) ! Compute zq, the height of Q relative to the polar stereographic plane: zq = erad * (cos(alpha) - 1.) ! Compute the parameter t which is the the distance ratio AQ:AR t = (erad2 + zq) / erad2 ! Compute xq and yq, the x and y coordinates of Q in polar stereographic space: xq = t * x yq = t * y ! Transform location of Q from (x,y,z) coordinates to (x3,y3,z3): x3q = x3p-xq*sinplon-yq*cosplon*sinplat+zq*cosplat*cosplon y3q = y3p+xq*cosplon-yq*sinplon*sinplat+zq*cosplat*sinplon z3q = z3p+yq*cosplat+zq*sinplat ! Compute the latitude and longitude of Q: qlon = atan2(y3q,x3q) / pi180 r3q = sqrt(x3q ** 2 + y3q ** 2) qlat = atan2(z3q,r3q) / pi180 return end subroutine xy_ll subroutine uevetouv(u,v,ue,ve,qlat,qlon,polelat,polelon) call ll_xy(qlat,qlon,polelat,polelon,x0,y0) call ll_xy(qlat,qlon+.1,polelat,polelon,x1,y1) angle = atan2(y1-y0,x1-x0) u = ue * cos(angle) - ve * sin(angle) v = ue * sin(angle) + ve * cos(angle) return end subroutine uevetouv subroutine uvtoueve(u,v,ue,ve,qlat,qlon,polelat,polelon) call ll_xy(qlat,qlon,polelat,polelon,x0,y0) call ll_xy(qlat,qlon+.1,polelat,polelon,x1,y1) angle = -atan2(y1-y0,x1-x0) ue = u * cos(angle) - v * sin(angle) ve = u * sin(angle) + v * cos(angle) return end subroutine uvtoueve