#include "dims.h" ! ! am_02/02: include general wave (origin Kelvin wave) C SUBROUTINE DT C **** C **** Advances T by one time step C **** C **** where T = temperature (deg K) C **** use input_module,only: iuivi,step,wave use init_module,only : iter use bndry_module,only: tb,tb2,tba,tbw,bnd,bnd2,bnda,bndw,ci use cons_module,only : len1,len2,len3,imax,imaxh,imaxp2,imaxp4, | kmax,kmaxm1,kmaxp1,dift,kut,rmassinv,expz,tbound,expzmid, | expzmid_inv,gask,grav,dtsmooth,dtsmooth_div2,dtx2inv,shapiro, | freq_semidi,tsurplus,avo,boltz,p0 implicit none #include "params.h" #include "fcom.h" #include "buff.h" #include "crates_tdep.h" #include "index.h" #include "phys.h" #include "vscr.h" ! ! Local: ! ! Bug fix 8/20/99: expt,expt2,expta were implicitly typed ! as reals instead of explicitly typed as complex. ! real:: wper complex :: expt,expt2,expta,exptw integer :: i,ntjp2k,ntjp1k,ntnmk,ntjm1k,ntjm2k,ncpk,nps1k,nps2k,k, | nqk,nqdhk,ntk,nuk,nvk,nuik,nvik,nlxxk,nlxyk,nlyxk,nlyyk,nkmk, | nmsk,nktk,nwtek,nwtik,nwk,ntnpk,ntnmnk,nxk,n real rstep ! C **** C **** /EDDDY/ contains arrays of DIFKT and DIFKK C **** DIFKT = eddy thermal conductivity used in subroutine DT. C **** DIFKK = eddy diffusion used in subroutines COMP, MINOR C **** and MINOR2 C **** C **** First calculate T(n+1) using a finite difference C **** approximation to the thermodynamic equation. The leap- C **** frog scheme is used with vertical diffusion treated C **** implicitly to second order accuracy. This leads to a C **** tridiagonal scheme requiring boundary conditions at the C **** top and bottom of the domain as implied by the C **** differential equation itself. The horizontal advection C **** is treated to 4th order accuracy. At level k the C **** differential equation becomes; C **** C **** P(k,n)*T(k-1,n+1) + Q(k,n)*T(k,n+1) + R(k,n)*T(k+1,n+1) C **** C **** = S(k,n) (1) C **** C **** where C **** k = vertical grid index, range k = 3/2, K-1/2, 1 C **** n = time-step index C **** C **** Variables are assumed to have indices (k,n) unless C **** otherwise stated. C **** _ C **** P(k,n) = -(G(k-1/2,n) - F(k-1/2,n)) | C **** | C **** Q(k,n) = G(k-1/2,n)+G(k+1/2,n)+F(k-1/2,n)- | C **** | C **** F(k+1/2,n)+cp*exp(-s)* | C **** | C **** (1/(2*Dt)+ai/cp+w*R/(cp*M)) | (2) C **** | C **** R(k,n) = -(G(k+1/2,n) + F(k+1/2,n)) | C **** | C **** S(k,n) = cp*exp(-s)*(T(k,n-1)/(2*Dt)- | C **** | C **** V.del(T(k,n))+(Q+QDH+GWT+QJ+QM)/cp- | C **** ae*exp(-s) | C **** _| C **** G = g*KT/(p0*H*Ds**2) C **** F = g*eps/(p0*2*H*Ds) C **** g = acceleration due to gravity (cm/sec**2) C **** KT = kT + H**2*rho*cp*kE C **** eps = kE*H**3*rho*g/T C **** kT = vertical molecular diffusion of heat C **** (ergs/cm/deg/sec) C **** kE = Eddy thermal diffusion (1/sec) C **** p = pressure (dynes/cm**2) C **** p0 = reference pressure (dynes/cm**2) C **** H = scale height = R*T/(M*g) (cm) C **** s = vertical pressure coordinate = ln(p0/p) C **** Ds = grid interval of vertical coordinate, s. C **** cp = specific heat at constant pressure C **** (ergs/deg/gm) C **** M = mean molecular mass C **** t = time (sec) C **** Dt = time step (sec) C **** ai = coefficient of implicit portion of radiative C **** cooling (1/sec) C **** ae = explicit portion of radiative cooling C **** (ergs/sec/gm) C **** w = ds/dt (1/sec) C **** V.del(T) = total 3d advection of T (deg/sec) C **** Q = total heating rate from all sources C **** (ergs/sec/gm) C **** QDH = heating due to 4th order horizontal C **** diffusion C **** GWT = heating due to gravity graves C **** rho = density (gm/cm**3) C **** R = gas constant per mole (ergs/deg) C **** C **** Boundary conditions: C **** C **** The lower boundary of the model is at k = 1 C **** T is given at this level as the sum of a global mean C **** value plus perturbations due to semidiurnal, diurnal C **** and annual tides. At level k = 3/2, (1) becomes: C **** C **** (Q(3/2,n)-P(3/2,n))*T(3/2,n+1) + R(3/2,n)*T(5/2,n+1) C **** C **** = S(3/2,n) - 2*P(3/2,n)*TB (3) C **** C **** The upper boundary is at k = K C **** At this level: C **** C **** (KT*d/ds + eps)T = 0 C **** C **** or in finite difference form: C **** C **** (G(K,n)-F(K,n))*T(K-1/2,n+1) - C **** C **** (G(K,n)+F(K,n))*T(K+1/2,n+1) = 0 C **** C **** This implies that at level k = K-1/2, (1) becomes: C **** C **** P(K-1/2,n)*T(K-3/2,n+1) + QT(K-1/2)*T(K-1/2,n+1) C **** C **** = S(K-1/2,n) (4) C **** C **** where C **** C **** QT(K-1/2,n) = G(K-1,n) + F(K-1,n)) + cp*exp(-s)* C **** C **** (1/(2*Dt)+ai/cp+w*R/(cp*M)) C **** C **** P(K-1/2,n) and S(K-1/2,n) are unchanged C **** C **** Equations (1) through (4) provide a tridiagonal system C **** of order (K-1), which is readily solved. C **** C **** C **** Evaluate tidal lower boundary condition for T. C **** This includes semidiurnal, diurnal and annual tides. C **** C **** T1 = temperature at lower boundary C **** C **** EXPT = exp(i*nusd*t) (time varying factor in C **** semidiurnal tide) C **** where i = sqrt(-1) C **** nusd = semidiurnal frequency (radians/sec) C **** t = current model time (sec) C **** rstep = float(step) EXPT = CEXP(CI*freq_semidi*rstep*ITER) C **** C **** EXPT2 = exp(i*nud*t) (time varying factor in C **** diurnal tide) C **** where nud = diurnal frequency (radians/sec) C **** EXPT2 = CEXP(CI*.5*freq_semidi*rstep*ITER) C **** C **** EXPTW = exp(i*nud*t) (time varying factor for wave) C **** where nud = frequency (radians/sec) C **** wper = wave(2) ! period of wave (days) exptw = cexp(ci*.5*freq_semidi*rstep/wper*iter) C **** C **** EXPTA = 1. (time variation factor for annual tide is 1. C **** since it is included in the latitudinal coefficients, C **** TBA, UBA and VBA) C **** EXPTA = 1. DO I = 1,LEN1 C **** C **** Semidiurnal tide C **** T1(I) = REAL(TB(J)*BND(I)*EXPT)+TBOUND C **** C **** where TB = latitudinal semdiurnal tidal variation C **** BND = longitudinal semdiurnal tidal variation C **** (TB and BND are calculated in BNDRY during C **** initialization) C **** EXPT = semidiurnal time variation C **** TBOUND = constant background temperature at C **** lower boundary C **** C **** Now add in the diurnal tide: C **** T1(I) = T1(I)+REAL(TB2(J)*BND2(I)*EXPT2) C **** C **** where TB2 = latitudinal diurnal tidal variation C **** BND2 = longitudinal diurnal tidal variation C **** (TB2 and BND2 are calculated in BNDRY2 during C **** initialization) C **** C **** C **** Now add in the general wave: C **** T1(I) = T1(I) + REAL(TBW(J)*BNDW(I)*EXPTW) C **** C **** where TBW = latitudinal wave variation C **** BNDW = longitudinal wave variation C **** (TBW and BNDW are calculated in BNDRY_WAVE during C **** initialization) C **** C **** Now add in the annual tide: C **** T1(I) = T1(I)+REAL(TBA(J)*BNDA(I)*EXPTA) C **** C **** where TBA = latitudinal variation C **** BNDA = longitudinal variation ( = 1. at present) C **** ENDDO ! write(6,"('dt: t1 (lower boundary, to be stored at kmaxp1):', ! | /(6e12.4))") t1 C **** C **** S1 = horizontal advection of T (temperature) (k+1/2) C **** C **** = (u/(r*cos(theta)))dT/d(lamda) + (v/r)dT/d(theta) C **** C **** (using 4th order difference approximation) C **** Note: The pointer NT, giving the position of the C **** temperature field in the buffer in the F-array, is C **** passed to the subroutine ADVEC. C **** Warning: ADVEC uses S11 through S15 as work space C **** CALL ADVEC(NT,S1) C **** C **** Add vertical advection of T to horizontal advection in S1 C **** C **** S1 = S1 + w*dT/ds (using 2nd order finite difference C **** approximation) (k+1/2) C **** C **** Here s = dimensionless vertical pressure coordinate C **** w = ds/dt C **** Warning: ADVECV uses S3 as work space C **** CALL ADVECV(NT,S1,T1) C **** C **** Shapiro smoother operating on T(n-1) in 2 stages C **** C **** 1) S2 = Result after meridional smoothing of T(n-1) C **** (k+1/2) C **** NTJP2K = NJP2+NTNM NTJP1K = NJP1+NTNM NTNMK = NJ+NTNM NTJM1K = NJM1+NTNM NTJM2K = NJM2+NTNM DO I = 1,LEN2 S2(I,1) = F(I,NTNMK)-shapiro*(F(I,NTJP2K)+F(I,NTJM2K)-4.* 1 (F(I,NTJP1K)+F(I,NTJM1K))+6.*F(I,NTNMK)) ENDDO C **** C **** 2) S3 = Result after zonal smoothing of S2 (k+1/2) C **** = smoothed T(n-1) C **** DO I = 3,LEN2-2 S3(I,1) = S2(I,1)-shapiro*(S2(I+2,1)+S2(I-2,1)-4.*(S2(I+1,1)+ | S2(I-1,1))+6.*S2(I,1)) ENDDO C **** C **** S1 = cp*exp(-s)*(V.del(T(n))-1/(2*DT))*T(n-1) (k+1/2) C **** C **** S8 = MBAR = 1./(psi(o2)/m(o2) + psi(o)/m(o) + C **** psi(n2)/m(n2)) (K+1/2) C **** NCPK = NCP-1 NPS1K = NJ + NPS NPS2K = NJ + NPS2 DO K = 1,KMAX NCPK = NCPK+1 NPS1K = NPS1K + 1 NPS2K = NPS2K + 1 DO I = 1,LEN1 S1(I,K) = .5*(F(I,NCPK)+F(I,NCPK+1))*expz(K)*(S1(I,K)-dtx2inv* 1 S3(I,K)) S8(I,K) = 1./(F(I,NPS1K)*rmassinv(1) + F(I,NPS2K)*rmassinv(2)+ 1 (1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(3)) ENDDO ENDDO C **** C **** Calculate total heat sources in S5 (k+1/2) C **** NQK = NQ NQDHK = NQDH DO I = 1,LEN2 C **** C **** First solar heating: C **** C **** S5(k) = .5*(Q(k)+Q(k+1)) (k+1/2) C **** S5(I,1) = .5*(F(I,NQK)+F(I,NQK+1)) C **** C **** Add heating from 4th order horizontal diffusion C **** (k+1/2) S5(I,1) = S5(I,1)+F(I,NQDHK) C **** ENDDO NTK = NJ + NT -1 NPS2K = NJ + NPS2 -1 DO K = 1,KMAX NTK = NTK + 1 NPS2K = NPS2K + 1 DO I = 1,LEN1 C **** (k+1/2) C **** Add heating due to atomic oxygen recombination C **** (k+1/2) C **** Reaction is: C **** O + O + M -> O2 + M C **** where M is N2, O2 or O C **** C **** Recombination rate = rk*n(O)**2*n (events/sec/cc) C **** C **** where rk = reaction rate C **** n(O) = number density of O (1/cc) C **** n = total number density (1/cc) C **** C **** Heating rate = h(O)*rk*n(O)**2*n/rho (ergs/sec/gm) C **** C **** where h(O) = heat generated per event (ergs) C **** rho = density (gm/cc) C **** In terms of mass mixing ratios, C **** C **** heating rate = h(O)*rk*((n*M)*psi(O)/M(O))**2*N0/M C **** C **** where M = mean molecular weight C **** M(O) = molecular weight of O C **** psi(O) = mass mixing of O C **** N0 = Avogadro's number C **** S5(I,K) = S5(I,K)+tsurplus*RKM12(I,K)*(p0*expz(K)*S8(I,K)/ 1 (boltz*F(I,NTK))*F(I,NPS2K)*rmassinv(2))**2*avo/S8(I,K) C **** ENDDO ENDDO C **** C **** Add Joule heating, QJ C **** C **** QJ = lamxx*(ui-u)**2 + (lamxy-lamyx)*(ui-u)* C **** (vi-v) + lamyy*(vi-v)**2 (k+1/2) C **** where: C **** lam-- are elements of ion drag tensor (1/sec) C **** u = zonal neutral velocity (cm/sec) C **** ui = zonal ion velocity (cm/sec) C **** v = meridional neutral velocity (cm/sec) C **** ui = meridional ion velocity (cm/sec) C **** C **** if(IUIVI = 2)then C **** S2 = -u (k+1/2) C **** else C **** S2 = ui -u C **** endif C **** S3 is set to (vi-v) or (-v) in the same way C **** NUK = NJ + NU NVK = NJ + NV NUIK = NUI NVIK = NVI IF(IUIVI.EQ.2)THEN DO I = 1,LEN2 S2(I,1) = -F(I,NUK) S3(I,1) = -F(I,NVK) ENDDO ELSE DO I = 1,LEN2 S2(I,1) = .5*(F(I,NUIK)+F(I,NUIK+1))-F(I,NUK) S3(I,1) = .5*(F(I,NVIK)+F(I,NVIK+1))-F(I,NVK) ENDDO ENDIF NLXXK = NLXX NLXYK = NLXY NLYXK = NLYX NLYYK = NLYY DO I = 1,LEN2 C **** C **** S4 = total Joule heating (k+1/2) C **** S4(I,1) = .5*(S2(I,1)**2*(F(I,NLXXK)+F(I,NLXXK+1))+S2(I,1)* 1 S3(I,1)*(F(I,NLXYK)+F(I,NLXYK+1)-F(I,NLYXK)-F(I,NLYXK+1))+ 2 S3(I,1)**2*(F(I,NLYYK)+F(I,NLYYK+1))) C **** C **** Add Joule heating to S5 (k+1/2) C **** S5(I,1) = S5(I,1)+S4(I,1)*1.5 C S5(I,1) = S5(I,1)+S4(I,1)*2. C S5(I,1) = S5(I,1)+S4(I,1)*3. C S5(I,1) = S5(I,1)+S4(I,1)*1. ENDDO c c Add solar heating QSOLAR and joule heating QJOULE to secondary histories: ! subroutine addfsech(name,long_name,units,f2d,londim,levdim, ! | nlev,lat) ! ! do k=1,kmaxp1 ! s10(imax+3,k) = s10(3,k) ! enddo ! call addfsech('QSOLAR','SOLAR HEATING from DT','DEG K', ! | s10,zimxp,zkmxp,zkmx,j) ! ! do k=1,kmaxp1 ! s4(imax+3,k) = s4(3,k) ! enddo ! call addfsech('QJOULE','JOULE HEATING from DT','DEG K', ! | s4,zimxp,zkmxp,zkmx,j) ! C **** C **** Add heating due to molecular diffusion, QM C **** C **** QM = Km/(rho*H**2)*((du/ds)**2+(dv/ds)**2) (ergs/gm) C **** C **** where Km = molecular viscosity (gm/cm/sec) (5) C **** C **** S10 = du/ds, S11 = dv/ds (k+1/2) C **** C **** Levels 2 (k = 5/2) through KMAX-1 (k = K-3/2) C **** C **** du/ds(k) = (u(k+1)-u(k-1))/(2*Ds) C **** dv/ds(k) = (v(k+1)-v(k-1))/(2*Ds) C **** NUK = NJ+NU NVK = NJ+NV DO K = 2,KMAXM1 NUK = NUK+1 NVK = NVK+1 DO I = 1,IMAXP4 S10(I,K) = (F(I,NUK+1)-F(I,NUK-1))/(2.*dz) S11(I,K) = (F(I,NVK+1)-F(I,NVK-1))/(2.*dz) ENDDO ENDDO C **** C **** Levels 1 (k = 3/2) and KMAX (k = K-1/2) C **** C **** At k = 3/2, we fit parabolas to the known values C **** of u and v at levels k = 1, k = 3/2 and k = 5/2. C **** We take the derivatives of these functions at C **** k = 3/2. C **** (Recall that level KMAXP1 contains the values of C **** u and v at level k = 1) C **** C **** At k = K-1/2, we make use of our upper boundary C **** condition (du/ds = dv/ds = 0 at level K) and fit a C **** parabola to levels (K+1/2), (K-1/2) and (K-3/2). C **** Taking u(K+1/2) = u(K-1/2) for the zero derivative C **** at level K, it can be shown that: C **** C **** du/ds(K-1/2) = 1/3*du/ds(K-3/2), C **** C **** and similarly for v. C **** NUK = NJ+NU NVK = NJ+NV DO I = 1,IMAXP4 S10(I,1) = (F(I,NUK) + 1./3.*F(I,NUK+1) - 4./3.*F(I,NUK+KMAX))/ 1 dz S11(I,1) = (F(I,NVK) + 1./3.*F(I,NVK+1) - 4./3.*F(I,NVK+KMAX))/ 1 dz C S10(I,1) = (F(I,NUK)-T2(I))/dz C S11(I,1) = (F(I,NVK)-T3(I))/dz c S10(I,1) = (F(I,NUK)+F(I,NUK+1)-2.*T2(I))/(2.*dz) C S11(I,1) = (F(I,NVK)+F(I,NVK+1)-2.*T3(I))/(2.*dz) S10(I,KMAX) = S10(I,KMAXM1)/3. S11(I,KMAX) = S11(I,KMAXM1)/3. ENDDO C **** C **** From (5), we have C **** C **** S10 = QM = g**2*M*Km/(p0*R*T*exp(-s))*(S10**2+S11**2) C **** (k+1/2) C **** NKMK = NKM-1 NTK = NJ+NT-1 DO K = 1,KMAX NKMK = NKMK+1 NTK = NTK+1 DO I = 1,IMAXP4 S10(I,K) = grav**2*S8(I,K)* 1 .5*(F(I,NKMK)+F(I,NKMK+1))/(p0*gask*expz(K)* 2 F(I,NTK))*(S10(I,K)**2+S11(I,K)**2) ENDDO ENDDO C **** C **** Add QM to other heat sources in S5 C **** DO I = 1,LEN2 S5(I,1) = S10(I,1)+S5(I,1) ENDDO c call addfsech('DT_S5',s5,zimxp,kmaxp1,kmaxp1,j) C **** C **** S1 = -cp*exp(-s)*(T(k,n-1)/(2*Dt) - V.del(T(k,n)) C **** C **** +Q/cp) C **** C **** = S1 - exp(-s)*S5 (k+1/2) C **** DO K = 1,KMAX DO I = 1,LEN1 S1(I,K) = S1(I,K) - expz(K)*S5(I,K) ENDDO ENDDO C **** C **** S2 = G = g*(kT + H**2*rho*cp*kE)/(p0*H*Ds**2) (k) C **** S3 = F = g*(kE*H**3*rho*g/T)/(p0*2*H*Ds) (k) C **** C **** S4 = H = R*T/(M*g) (k) C **** S5 = rho = p0*exp(-s)*M/(R*T) (k) C **** S6 = T (k) C **** C **** Levels 2 through (K-1) C **** NTK = NJ+NT-1 NMSK = NJ+NMS-1 DO K = 2,KMAX NTK = NTK+1 NMSK = NMSK+1 DO I = 1,LEN1 S6(I,K) = .5*(F(I,NTK)+F(I,NTK+1)) S4(I,K) = gask*S6(I,K)/F(I,NMSK+1) S5(I,K) = p0*expzmid_inv*expz(K)/S4(I,K) S4(I,K) = S4(I,K)/grav ENDDO ENDDO C **** C **** Levels 1 and K C **** C **** (recall TB in T(KMAXP1), T(K) = T(K-1/2) ) C **** NTK = NJ+NT NMSK = NJ+NMS DO I = 1,LEN1 S6(I,1) = F(I,NTK+KMAX) S6(I,KMAXP1) = F(I,NTK+KMAXM1) S4(I,1) = gask*S6(I,1)/F(I,NMSK) S4(I,KMAXP1) = gask*S6(I,KMAXP1)/F(I,NMSK+KMAX) S5(I,1) = p0*expzmid_inv*expz(1)/S4(I,1) S5(I,KMAXP1) = p0*expzmid*expz(KMAX)/S4(I,KMAXP1) S4(I,1) = S4(I,1)/grav S4(I,KMAXP1) = S4(I,KMAXP1)/grav ENDDO C **** C **** S2 = G (k) C **** S3 = F (k) C **** NKTK = NKT-1 NCPK = NCP-1 DO K = 1,KMAX NKTK = NKTK+1 NCPK = NCPK+1 DO I = 1,LEN1 S2(I,K) = grav*(F(I,NKTK)+S4(I,K)**2*S5(I,K)*F(I,NCPK)* 1 DIFT(K))/(p0*S4(I,K)*dz**2) S3(I,K) = grav**2*DIFT(K)*S4(I,K)**2*S5(I,K)/ 1 (S6(I,K)*p0*2.*dz) ENDDO ENDDO C **** C **** S9 = -P (k+1/2) C **** S10 = -Q (k+1/2) C **** S11 = -R (k+1/2) C **** S12 = -S (k+1/2) C **** C **** (See (2) above for definitions of P, Q, R and S) C **** C **** Levels 3/2 through K-3/2 C **** DO I = 1,LEN2 - LEN1 S9(I,1) = S2(I,1) - S3(I,1) S10(I,1) = -S2(I,1) - S2(I,2) - S3(I,1) + S3(I,2) S11(I,1) = S2(I,2) + S3(I,2) S12(I,1) = S1(I,1) ENDDO C **** C **** At level K-1/2, we use (4) and have C **** DO I = 1,LEN1 S9(I,KMAX) = S2(I,KMAX) - S3(I,KMAX) S10(I,KMAX) = -S2(I,KMAX) - S3(I,KMAX) S11(I,KMAX) = 0. S12(I,KMAX) = S1(I,KMAX) ENDDO C **** C **** C **** C **** S10 = S10 - cp*exp(-s)*(1/(2*DT)+ai/cp+w*R/(cp*M)) C **** to complete -Q (k+1/2) C **** C **** S12 = S12 + exp(-s)*ae (k+1/2) C **** to complete -S C **** NWTEK = NWTE NWTIK = NWTI NCPK = NCP NWK = NJNP+NW DO I = 1,LEN2 C **** C **** S13 = cp*(1/(2*Dt)+ai/cp+w*R/(cp*M)) (k+1/2) C **** S13(I,1) = .5*(F(I,NCPK)+F(I,NCPK+1))*(dtx2inv+F(I,NWTIK))+ 1 .5*(F(I,NWK)+F(I,NWK+1))*gask/S8(I,1) S12(I,1) = S12(I,1)+F(I,NWTEK) ENDDO DO K = 1,KMAX DO I = 1,LEN1 C **** C **** S10 = S10 - exp(-s)*S13 C **** S10(I,K) = S10(I,K)-expz(K)*S13(I,K) ENDDO ENDDO C **** C **** Set boundary conditions C **** DO I = 1,LEN1 C **** C **** Lower boundary at k=3/2 C **** S10 = -(Q - P) = S10 - S9 C **** S12 = -(S-2*P*TB) = S12 - 2*S9*T1 C **** S9 = P = 0. C **** S10(I,1) = S10(I,1) - S9(I,1) S12(I,1) = S12(I,1) - 2.*S9(I,1)*T1(I) S9(I,1) = 0. ENDDO C **** C **** Call TRSLOV to solve tridiagonal system for T(n+1). C **** Solution is placed in buffer NJNP C **** Note: S13, S14 and S15 are used as work space although C **** only S15 is passed as a parameter. Use is made of the C **** fact that they are contiguous in /VSCR/. C **** NTNPK = NJNP+NT CALL TRSOLV(S9,S10,S11,S12,S15,F(1,NTNPK),LEN1,3,LEN1-2,KMAXP1, 1 1,KMAX,1) C **** C **** C **** Filter T C **** IF (KUT(J).LT.IMAXH) THEN NTNPK = NJNP+NT CALL FILTER(NTNPK,KMAX,KUT(J)) ENDIF C **** C **** Time smoothing of T C **** NTNMK = NJ+NTNM NTK = NJ+NT NTNPK = NJNP+NT NTNMNK = NJNP+NTNM DO I = 1,LEN2 C **** C **** T(smooth) = TNMN = .5*(1.-alpha)*(TNM+TNP) + alpha*TN C **** F(I,NTNMNK) = dtsmooth_div2*(F(I,NTNMK)+F(I,NTNPK)) + | dtsmooth*F(I,NTK) ENDDO C **** C **** Store lower boundary values in level KMAXP1 C **** NTNPK=NJNP+NT+KMAX DO I=1,IMAXP4 F(I,NTNPK)=T1(I) ENDDO C **** C **** Insert periodic points for updated T. C **** NXK=NJNP+NT-1 DO N = 1,2 DO I = 1,2 DO K=1,KMAXP1 F(I,NXK+K)=F(IMAX+I,NXK+K) F(IMAXP2+I,NXK+K)=F(2+I,NXK+K) ENDDO ENDDO NXK=NJNP+NTNM-1 ENDDO C **** C **** SET MINIMUM VALUE OF TN AT 100K C **** NTK = NJNP+NT DO I=1,LEN3 if (f(i,ntk) < 100.) f(i,ntk) = 100. ENDDO RETURN END C