#include "dims.h" SUBROUTINE NMCBND(istep,mtime) use cons_module,only: len1,imaxp4,jmax,atm_amu,grav,gask,dt use init_module,only: iyear,iday,iter implicit none ! integer,intent(in) :: istep,mtime(4) C **** C **** Define Z (and T) lower boundary from NMC data: C **** (called from early in advnce after day advance) C **** #include "params.h" #include "ncep.h" real,save :: glon(zimxp),d4(zimxp,zjmx,npr,4) real :: SDIF(3) integer,save :: iyd4(4) real :: PRESS(3) = (/30.,10.,5./) ! c data mssfile/'/FOSTER/nmc/79001-93120.dat '/ c data iyd0,iyd1/79001,93120/ ! first and last days of nmc data ! data mssfile/'/FOSTER/nmc/79001-94151.dat '/ ! updated 3/9/95 ! data iyd0,iyd1/79001,94151/ ! updated 3/9/95 ! ! 6/9/99: mssfile 79001-98120.dat was made from ieee formatted files from ! /WUF/nmc (Fei Wu). See ~foster/timegcm/nmc. ! ! mssfile updated 6/9/99 character(len=32),parameter :: | mssfile = '/FOSTER/nmc/79001-98120.dat ' ! updated 6/9/99 integer :: iyd0=79001, iyd1=98120 ! updated 6/9/99 ! ! character(len=32),parameter :: ! | mssfile = '/FOSTER/nmc/98116-98120.dat ' ! test file ! integer :: iyd0=98116, iyd1=98120 ! for test file ! integer :: mxrec=7060, lu=98, idatsec=0, iydp=99999 integer :: i,iyd,iutsec,iprint,ih,im,is,k,j,n,ipos,ier integer,external :: inciyd ! To save z and t on secondary histories (redundant in pressure): real :: fmin,fmax,znmcik(zimxp,zkmxp),tnmcik(zimxp,zkmxp) ! ! data lu/98/, idatsec/0/, iydp/99999/ c c Determine current calendar day and time in hrs,min,sec: c C write(6,"('enter nmcbnd: iter=',i6,' istep=',i6)") iter,istep iyd = (iyear-1900)*1000+iday if (iyd.lt.iyd0.or.iyd.gt.iyd1-1) then write(6,"(' ')") write(6,"('>>> NMCBND: iyd outside NMC data ', + 'range: desired iyd=',i5,' NMC data available for ',i5, + ' to ',i5)") iyd,iyd0,iyd1 stop 'NMCBND' endif iutsec = ifix(amod(float(iter)*dt,86400.)) call isec2hms(iutsec,ih,im,is) c c If first iteration, get 4 days of nmc data from which time c interpolations start (these will be advanced if model c crosses a day boundary): c Also put message out to stdout. Let nmczday do sample output c to stdout if iprint==1 (1st day only) c ! iprint = 0 iprint = 1 if (istep==1) then write(6,"(/72('-'))") write(6,"('NMCBND:')") write(6,"('Lower boundary for heights (10mb) will be obtained ', + ' at each time step from')") write(6,"(' NMC data on mss unformatted direct access file ', + a)") mssfile(1:len_trim(mssfile)) write(6,"('Interpolation starts at yyddd ',i5,' hr:min:sec ', + i2,':',i2,':',i2,' (iter=',i8,')'/)") iyd,ih,im,is,iter do i=1,imaxp4 glon(i) = -190.+(i-1)*5. enddo do i=1,4 ! iprint = 0 iprint = 1 ! if (i == 1) iprint = 1 iyd4(i) = inciyd(iyd,i-2) call nmczday(lu,mssfile,mxrec,iyd0,iyd1,iyd4(i),d4(1,1,1,i), | glon,imaxp4,jmax,npr,ier,iprint) if (ier.ne.0) then write(6,"(' ')") write(6,"('>>> ERROR FROM NMCZDAY (called from NMCBND', + ' at first iteration)')") write(6,"(' ')") stop 'NMCZDAY' endif enddo write(6,"(72('-')/)") else c c Advance the 4-day interpolates if necessary: c if (iyd.gt.iydp) then do i=1,4 iyd4(i) = inciyd(iyd4(i),iyd-iydp) ! presumably iyd-iydp=1 if (i.le.3) d4(:,:,:,i) = d4(:,:,:,i+1) enddo call nmczday(lu,mssfile,mxrec,iyd0,iyd1,iyd4(4),d4(1,1,1,4), | glon,imaxp4,jmax,npr,ier,iprint) if (ier.ne.0) then write(6,"(' ')") write(6,"('>>> ERROR FROM NMCZDAY (called from NMCBND)', + /' was advancing day from ',i5,' to ',i5)") iydp,iyd write(6,"(' ')") stop 'NMCZDAY' endif endif endif c c Do interpolation to current time, defining znmc: c ipos = 2 if (iyd.eq.iyd0) ipos = 1 if (iyd.eq.iyd1-1) ipos = 3 do k=1,npr do j=1,jmax call timeterp(d4(1,j,k,1),d4(1,j,k,2),d4(1,j,k,3),d4(1,j,k,4), + imaxp4,ipos,idatsec,iutsec,znmc(1,j,k),0,ier) enddo enddo C **** C **** Calculate T at lower boundary from: C **** C **** T = MBAR*g/R * dZ/ds C **** C **** and C **** C **** dZ/ds = ((ln(p2/p3))**2*(z2-z1) + (ln(p1/p2))**2* C **** (z3-z1))/(-ln(p1/p2)*ln(p21/p32)*ln(p31/p12)) C **** C **** Calculate differences between the three pressure levels C **** in terms of the pressure coordinate, s. C **** DO N = 1,3 SDIF(N) = ALOG(PRESS(N)/PRESS(1+MOD(N,3))) ENDDO C **** C **** Calculate C **** C **** DO J = 1,JMAX DO I = 1,LEN1 TNMC(I,J) = atm_amu*grav/gask*(SDIF(2)**2* 1 (ZNMC(I,J,2)-ZNMC(I,J,1)) + SDIF(1)**2* 2 (ZNMC(I,J,3)-ZNMC(I,J,2))) / 3 (-SDIF(1)*SDIF(2)*SDIF(3)) ENDDO ENDDO ! ! Save output arrays zncep(zimxp,zjmx,npr),tncep(zimxp,zjmx) on ! secondary histories (redundant in pressure). ! ! call fminmax(znmc(:,:,2),zimxp*zjmx,fmin,fmax) ! write(6,"('istep=',i5,' mtime=',3i4,' nmc Z at 10 mb min,max=', ! | 2e12.4)") istep,mtime(1:3),fmin,fmax ! call fminmax(tnmc,zimxp*zjmx,fmin,fmax) ! write(6,"(' nmc T at 10 mb min,max=',2e12.4)") fmin,fmax ! do j=1,zjmx ! do k=1,zkmxp ! redundant in pressure for secondary histories ! znmcik(:,k) = znmc(:,j,2) ! save at 10 mb level ! enddo ! call addfsech('Z_NMC',' ',' ',znmcik,zimxp,zkmxp,zkmxp,j) ! ! do k=1,zkmxp ! redundant in pressure for secondary histories ! tnmcik(:,k) = tnmc(:,j) ! tnmc is at 10 mb ! enddo ! call addfsech('T_NMC',' ',' ',tnmcik,zimxp,zkmxp,zkmxp,j) ! enddo ! iydp = iyd return end