#include "dims.h" SUBROUTINE DT use input_module,only: iuivi,ncep,nmc,step use init_module,only: iter use bndry_module,only: tb,tb2,tba,bnd,bnd2,bnda,ci use cons_module,only: kut,len1,len2,len3,kmax,kmaxp1,kmaxm1,expz, | imax,imaxh,imaxp4,imaxp2,freq_semidi,freq_3m3,shapiro,dtx2inv, | grav,p0,gask,expzmid,expzmid_inv,dtsmooth,dtsmooth_div2 implicit none C **** C **** Advances T by one time step C **** C **** where T = temperature (deg K) C **** #include "params.h" #include "fcom.h" #include "buff.h" #include "cons3m3.h" #include "index.h" #include "phys.h" #include "zatmos.h" #include "vscr.h" #include "cool.h" #include "ncep.h" 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 **** #include "diffk.h" ! ! Local: real :: RADCOOL(ZIMXP,ZKMXP),QSOLAR(ZIMXP,ZKMXP),rstep integer :: n,i,k,ntjp2k,ntjp1k,ntnmk,ntjm1k,ntjm2k,ncpk,nqk, | nmshk,nmsk,nnmbark,nps1k,nps2k,nlxxk,nlyyk,nlxyk,nlyxk, | nuik,nvik,nuk,nvk,nqdhk,ngwvtk,nkmk,nktk,ntk,nwtek,nwtik, | nwk,ntnpk,nxk,ntnmnk C **** C **** Complex variables for the time phase of each tidal C **** contribution C **** COMPLEX EXPT, EXPT2, EXPTA, EXPT3M3 C **** C **** Where: C **** C **** EXPT = exp(i*nu*dt*iter) C **** nu = semidiurnal frequency (radians/sec) C **** dt = timestep (seconds) C **** iter = time-step count since time zero C **** similarly, C **** C **** EXPT2 refers to the diurnal frequency C **** C **** EXPTA is the time factor for the annual tide C **** C **** EXPTM3M is the quantity needed for the (3,-3) C **** for the 2-day planetary wave C **** C **** 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+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 **** C **** Evaluate tidal lower boundary condition for T. C **** This includes semidiurnal, diurnal, and annual tides C **** plus '2-day' wave 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 **** 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. C **** C **** EXPT3M3 = exp(i*nu3m3*t) (time varying factor for C **** (3,-3) wave) C **** where nu3m3 = frequency of '2-day' wave (radians/sec) C **** EXPT3M3 = CEXP(CI*freq_3m3*rstep*ITER) C **** C TN boundary: C **** Semidiurnal tide C T1(I) = REAL(TB(J)*BND(I)*EXPT)+TBOUND C T1(I) = REAL(TB(J)*BND(I)*EXPT) C **** 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 **** 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 **** Now add in the annual tide: C **** C T1(I) = T1(I)+REAL(TBA(J)*BNDA(I)*EXPTA) C T1(I) = T1(I)+ZATMOS_TN(J) C **** C **** where TBA = latitudinal variation C **** BNDA = longitudinal variation ( = 1. at present) C **** C **** Add in (3,-3) wave C **** C T1(I) = T1(I)+REAL(TB3M3(J)*BND3M3(I)*EXPT3M3) C **** C **** where TB3M3 = latitudinal variation of (3,-3) wave C **** BND3M3 = longitudinal variation of (3,-3) wave C **** (TB3M3 and BND3M3 are calculated in BNDRY3M3 C **** during initialization) C **** ! ! If ncep==1, then tncep is ncep boundary, otherwise (ncep==0), use ! boundary from zatmos_tn: ! if (ncep <= 0 .and. nmc <= 0) THEN DO I = 1,LEN1 T1(I) = REAL(TB(J)*BND(I)*EXPT) T1(I) = T1(I)+REAL(TB2(J)*BND2(I)*EXPT2) T1(I) = T1(I)+ZATMOS_TN(J) T1(I) = T1(I)+REAL(TB3M3(J)*BND3M3(I)*EXPT3M3) ENDDO else if (ncep > 0) then DO I = 1,LEN1 T1(I) = REAL(TB(J)*BND(I)*EXPT)+tncep(i,j) ! new ncep T1(I) = T1(I)+REAL(TB2(J)*BND2(I)*EXPT2) T1(I) = T1(I)+REAL(TB3M3(J)*BND3M3(I)*EXPT3M3) ENDDO elseif (nmc > 0) then DO I = 1,LEN1 T1(I) = REAL(TB(J)*BND(I)*EXPT)+tnmc(i,j) ! old nmc T1(I) = T1(I)+REAL(TB2(J)*BND2(I)*EXPT2) T1(I) = T1(I)+REAL(TB3M3(J)*BND3M3(I)*EXPT3M3) ENDDO endif endif 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 **** NCPK = NCP-1 DO K = 1,KMAX NCPK = NCPK+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)) ENDDO ENDDO C **** C **** Calculate total heat sources in S5 (k+1/2) C **** NQK = NQ NMSHK =NMSH NMSK = NJ+NMS NNMBARK = NNMBAR NPS1K = NJ+NPS NPS2K = NJ+NPS2 NLXXK = NLXX NLYYK = NLYY NLXYK = NLXY NLYXK = NLYX NUIK = NUI NVIK = NVI NUK = NJ+NU NVK = NJ+NV NQDHK = NQDH NGWVTK = NGWVT 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 Save solar heating in s10 for use below: c s10(i,1) = s5(i,1) C **** C **** Add heating from 4th order horizontal diffusion C **** (k+1/2) S5(I,1) = S5(I,1)+F(I,NQDHK) C **** C **** Add heating due to gravity waves C **** S5(I,1) = S5(I,1)+F(I,NGWVTK) 53 CONTINUE 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 **** 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 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 **** C S5(I,1) = S5(I,1)+S4(I,1)*3. S5(I,1) = S5(I,1)+S4(I,1)*2. 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: c subroutine addfsech(fname,f,idim1,idim2,ilat) c ! ncpk = ncp ! do i=1,len2 ! qsolar(i,1) = s10(i,1)*86400./f(i,ncpk) ! enddo ! call addfsech('QSOLAR',' ',' ',qsolar,zimxp,kmaxp1,kmax,j) ! ! call addfsech('F_NUI' ,' ',' ',f(1,nui) ,zimxp,kmaxp1,kmax,j) ! call addfsech('F_NVI' ,' ',' ',f(1,nvi) ,zimxp,kmaxp1,kmax,j) ! call addfsech('F_NU' ,' ',' ',f(1,nj+nu),zimxp,kmaxp1,kmax,j) ! call addfsech('F_NV' ,' ',' ',f(1,nj+nv),zimxp,kmaxp1,kmax,j) ! call addfsech('F_NLXX',' ',' ',f(1,nlxx) ,zimxp,kmaxp1,kmax,j) ! call addfsech('F_NLXY',' ',' ',f(1,nlxy) ,zimxp,kmaxp1,kmax,j) ! call addfsech('F_NLYX',' ',' ',f(1,nlyx) ,zimxp,kmaxp1,kmax,j) ! call addfsech('F_NLYY',' ',' ',f(1,nlyy) ,zimxp,kmaxp1,kmax,j) ! call addfsech('DT_S4' ,' ',' ',s4 ,zimxp,kmaxp1,kmax,j) ! **** 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 **** NMSHK = NMSH-1 NKMK = NKM-1 NTK = NJ+NT-1 DO K = 1,KMAX NMSHK = NMSHK+1 NKMK = NKMK+1 NTK = NTK+1 DO I = 1,IMAXP4 S10(I,K) = grav**2*F(I,NMSHK)* 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 **** 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 DIFKT(I,K))/(p0*S4(I,K)*dz**2) S3(I,K) = grav**2*DIFKT(I,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+w*R/(cp*M)) (k+1/2) C **** to complete -Q C **** C **** S12 = S12 + exp(-s)*ae (k+1/2) C **** to complete -S C **** NWTEK = NWTE NWTIK = NWTI NCPK = NCP NMSHK = NMSH 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/F(I,NMSHK) 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 addfsech('DT_S9' ,' ',' ',s9 ,zimxp,kmaxp1,kmax,j) ! call addfsech('DT_S10',' ',' ',s10,zimxp,kmaxp1,kmax,j) ! call addfsech('DT_S11',' ',' ',s11,zimxp,kmaxp1,kmax,j) ! call addfsech('DT_S12',' ',' ',s12,zimxp,kmaxp1,kmax,j) ! call addfsech('DT_NJTN',' ',' ',f(1,nj+nt),zimxp,kmaxp1,kmax,j) CALL TRSOLV(S9,S10,S11,S12,S15,F(1,NTNPK),LEN1,3,LEN1-2,KMAXP1, 1 1,KMAX,1) C **** C **** Store lower boundary values in level KMAXP1 C **** NTNPK=NJNP+NT+KMAX DO I=1,IMAXP4 F(I,NTNPK)=T1(I) ENDDO ! call addfsech('DTF_NT1' ,' ',' ',f(1,njnp+nt),zimxp,kmaxp1,kmax,j) C **** C **** Filter T C **** IF (KUT(J).LT.IMAXH) THEN NTNPK = NJNP+NT CALL FILTER(NTNPK,KMAXP1,KUT(J)) ENDIF ! call addfsech('DTF_NT2' ,' ',' ',f(1,njnp+nt),zimxp,kmaxp1,kmax,j) C **** C **** Insert periodic points for updated T. ! 6/10/98: moved periodic points from bottom to here, so they ! can be referenced in time smoothing. 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 **** Time smoothing of T C **** NTNMK = NJ+NTNM NTK = NJ+NT NTNPK = NJNP+NT NTNMNK = NJNP+NTNM DO I = 1,LEN3 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 ! call addfsech('DTF_NT3' ,' ',' ',f(1,njnp+nt),zimxp,kmaxp1,kmax,j) ! ! tstadj and adjust appear to be problematic starting at 2,9,52 ! (Model runs through the problem successfully if this call is removed) CALL TSTADJ ! ! Save temperature: ! call addfsech('NTN',' ',' ',f(1,njnp+nt),zimxp,kmaxp1,kmax,j) ! call addfsech('DTF_NT4' ,' ',' ',f(1,njnp+nt),zimxp,kmaxp1,kmax,j) c c Save cooling terms (in deg K per day) c Total radiative cooling: c nwtek = nwte-1 nwtik = nwti-1 ncpk = ncp-1 ntnpk = njnp+nt-1 do k=1,kmax nwtek = nwtek+1 nwtik = nwtik+1 ncpk = ncpk+1 ntnpk = ntnpk+1 do i=1,imaxp4 s14(i,k) = f(i,nwtek) / expz(k) + f(i,nwtik) * + f(i,ncpk) * f(i,ntnpk) enddo enddo ncpk = ncp do i=1,len2 radcool(i,1) = s14(i,1) * 86400./f(i,ncpk) enddo ! call addfsech('RADCOOL',' ',' ',RADCOOL,zimxp,kmaxp1,kmax,j) RETURN END