SUBROUTINE STREAMF(im10) c c file streamf_base8sr.f FORTRAN routine TO COMPUTE 2D MODEL RESID CIRC FROM TEMPS, HEAT RATES, C by solving stream function as in the Garcia-Solomon model, and the 2D coupled model C solves by SUCCESSIVE OVER-RELAXATION, solves for 36 10-day means C C This also includes extra X (Chi) term not contained in GS83 formulation, and also small modifications C to the Czz, and Cz coefficients. HOWever, these changes made ONLY VERY SMALL DIFFERENCES in the C Chi, w, and v fields. Note also that the solution converged a little FASTER with these changes C C Note: using constant N2 has very little effect except in tropical troposphere where C circ cells are 2-3X bigger with n2 from temps (with old heating rates) C C ***** NOTE THE COPIOUS DOCUMENTATION **************** C C ALL CALCULTIONS DONE IN MKS UNITS ==> v*, w*, Kyy, Kzz converted to cm(cm2)/sec C *** USES dely and delz, and rho59, identical to that used in transport routines *** C *** INCLUDES routine KWAVE for SAO ****** C c *** EVERYTHING IS IN DOUBLE PRECISION (REAL*8) and FOR 59 levels *********** #include "com2d.h" INTEGER ITERC C if it doesn't converge after 500 iterations, it stops PARAMETER(ITERC=500) REAL*8 srr, rr1, rad, om1, hh, kap, xpi, qp1, dely, dely2 REAL*8 delz, delz2, pres59(59), zz59(59), rho59(59), gr59(59) REAL*8 sf1, sf2, sf3, sf4, sf6, sf7, sf8, sf9, sf0, sf5 REAL*8 adiffc, dtest, kk1, kk2 REAL*8 xlat5(37), sinl5(37), cosl5(37), tana(37), ff0(37) REAL*8 ff1(37,59), ff2(37,59), kru(37,59), kyyad, qvad REAL*8 kry(37,59), kry1(37,59), df1dy(37,59), kyysh REAL*8 heat1(37,59), heat1s(37,59), lh1(37,59), kyy(37,59) REAL*8 temp1(37,59), dtt1(37,59), dt1y(37,59), kyymeso(37,59) REAL*8 ddh(37,59), ddhs(37,59), bbc2(37), fxall(37,59), qv(37,59) REAL*8 qtot(37,59), qtots(37,59), qtots2(37,59), c0(37,59) REAL*8 czz(37,59), cz(37,59), czy(37,59), cy(37,59), cyy(37,59) REAL*8 cff(37,59), cffs(37,59), chi0(37,59,2), chif(37,59) REAL*8 dcdy(37,59), dcdz(37,59), wchi(37,59), diffc(37,59) REAL*8 dudy(37,59), dtdy(37,59), dtdz(37,59), lhex(37,59) REAL*8 dt1z(37,59), dt1yy(37,59), dhdy(37,59), dfxdz(37,59) REAL*8 XN2(37,59), UBAR1(37,59), DUDZ(37,59), XKZZT(37,59) REAL*8 FXTT(37,59), FK(37,59) REAL*8 w1(18,59), WBAR(2,59), W(18,59), V(19,58) COMMON/NEWSF1/XN2, UBAR1, DUDZ, XKZZT, FXTT, FK COMMON/NEWSFC/rad, om1, delz, hh, rr1, kap, xlat5, pres59, c zz59, rho59, gr59 c keep at REAL*4 -- maybe reduce these to 46 levels after testing is done..... COMMON/DYNEXTRA/ubarall(37,59,36), kzzall37(37,59,36), * fxtall(37,59,36), qvall(37,59,36), kruall(37,59,36), * fkall(37,59,36), walla(37,59,36), valla(37,59,36), * heatall(37,59,36), ddhall(37,59,36), lhall(37,59,36), * qtotall(37,59,36), cffall(37,59,36), chiall(37,59,36) c C srr=successive over-relaxation parameter: 1= stan (no over-relxtion) 1 Dont need to convert C heat10 in units of deg K/day => convert to K/sec C lh10 in units of deg K/day => convert to K/sec C Kyy10 in units of m2/sec = Don't need to convert C double check current limit of 5.e10 cm2/s (this can be reduced) C set minimum Kyy to 1.e8 cm2/s (1.e6 seems to small) C but note that only 1.e10 is the Kyy limit at the top levels c NOTE: Kyy and Kzz are used here only to compute ddh term c NMC-based input has been lowered by 1/2 grid point, and level 59 = level 58 c do 105 ik=1,59 do 105 ij=1,37 temp1(ij,ik) = temp10(ij,ik,im10) ubar1(ij,ik) = ubar10(ij,ik,im10) qv(ij,ik) = qv10(ij,ik,im10) c Make Delf>0.0 == 0.0 if (qv(ij,ik) .gt. 0.) qv(ij,ik) = 0.0 c if (qv(ij,ik) .gt. 0.) qv(ij,ik) = qv(ij,ik)*cosl5(ij) lh1(ij,ik) = 0.0 ddh(ij,ik) = 0.0 c initialize lh1 array here before defining below kyy(ij,ik) = kyy10(ij,ik,im10) if (kyy(ij,ik) .lt. 1.d4) kyy(ij,ik) = 1.d4 if (kyy(ij,ik) .gt. 5.d6) kyy(ij,ik) = 5.d6 105 continue c c adjust Kyy, Delf in tropical lower troposphere (levs 1-4) to get correct interhemispheric transport C linearly interpolate DELF between 20S-20N, set Kyy to a minimum of 150 (e8) cm2/sec C blend over level 4 to get a smooth transistion. Don't worry about consistency with QBARY, C since this is just the TROPOSPHERE! C kyyad = 150.d4 do 111 ik=1,4 do 111 ij=1,37 if (kyy(ij,ik) .lt. kyyad) then if (ik .le. 3) kyy(ij,ik) = kyyad if (ik .eq. 4) kyy(ij,ik) = (kyyad + kyy(ij,ik))/2. endif 111 CONTINUE do 112 ik=1,4 do 112 ij=16,22 qvad = (qv(23,ik)-qv(15,ik))/8.*(ij-15) + qv(15,ik) if (ik .le. 3) qv(ij,ik) = qvad if (ik .eq. 4) qv(ij,ik) = (qv(ij,ik) + qvad)/2. 112 continue c Joan's INPUT heating rates still at 58 levs , set level 59 = 58 do 206 ik=1,58 do 206 ij=1,37 206 heat1(ij,ik) = heat10(ij,ik,im10)/86400. do 107 ij=1,37 107 heat1(ij,59) = heat1(ij,58) do 116 ik=1,7 do 116 ij=1,37 116 lh1(ij,ik) = lh10(ij,ik,im10)/86400. c c add in extra latent heating mainly in sub-tropics, so Hadley cell is not as strong, and more spread out c NO, DON'T DO THIS WITH PPM c do 390 ik=1,59 c do 390 ij=1,37 c 390 lhex(ij,ik) = 0.0 c c do 391 ij=7,31 c 391 lhex(ij,3) = 2.0/86400.*dcosd(3.*(dabs(xlat5(ij))-30.)) c c do 392 ij=9,29 c 392 lhex(ij,4) = 2.0/86400.*dcosd(3.6*(dabs(xlat5(ij))-25.)) c c do 393 ij=11,27 c 393 lhex(ij,5) = 1.35/86400.*dcosd(4.5*(dabs(xlat5(ij))-20.)) c c do 394 ij=13,25 c 394 lhex(ij,6) = .65/86400.*dcosd(6.*(dabs(xlat5(ij))-15.)) cc add to latent heat c do 398 ik=1,59 c do 398 ij=1,37 c 398 lh1(ij,ik) = lh1(ij,ik) + lhex(ij,ik) c C background Rayleigh fric w/ latitude dependence only, 1/RF = ~100 days in tropics, C 50 days at mid-high lats, decrease to 1/300 days in troposphere everywhere since we already C have momentum source from qv from the NMC data - convert to 1/sec do 106 ik=1,59 do 106 ij=1,37 106 kry1(ij,ik) = 150.-100*dexp(-dcosd(xlat5(ij))**4) do 415 ik=1,7 do 415 ij=1,37 415 kry1(ij,ik) = 300. do 416 ij=1,37 kry1(ij,8) = .75*kry1(ij,7) + .25*kry1(ij,11) kry1(ij,9) = .5*kry1(ij,7) + .5*kry1(ij,11) 416 kry1(ij,10) = .25*kry1(ij,7) + .75*kry1(ij,11) do 419 ik=1,59 do 419 ij=1,37 419 kry(ij,ik) = 1./kry1(ij,ik)/86400. c C ; Now compute derivatives of temp1, ubar1, dudy, dudz, dtdy, dtdz = (37,59) call derv4(1, ubar1, dudy, dely, L$1, Z$59, 1) call derv4(1, ubar1, dudz, delz, L$1, Z$59, 2) call derv4(1, temp1, dtdy, dely, L$1, Z$59, 1) call derv4(1, temp1, dtdz, delz, L$1, Z$59, 2) C Brunt-Visala freq as defined in Lindzen (1981), with minimum value set to 3.2e-4 to keep strmfunc elliptic do 130 ik=1,59 do 130 ij=1,37 xn2(ij,ik) = (dtdz(ij,ik) + gr59(ik)/1004.)*gr59(ik)/temp1(ij,ik) if (xn2(ij,ik) .lt. 3.2d-4) xn2(ij,ik) = 3.2d-4 130 continue C call gravity wave routine to compute Kzzt (m2/sec) and Fxtt (m/sec2) C after dudz, xn2 have been determined CALL GWAVE(im10) c call Kelvin wave routine to get FK(37.59) (m/sec2) momentum forcing for SAO CALL KWAVE(im10) C C C NOW COMPUTE FORCING TERMS: C C First the Dh term (eq 7) in GS83, eddy diffusion of heat by GWs, keep sign convention as they have it C although the sign of ddh here makes little difference, and the term in general seems to have a C NEGLIGIBLE EFFECT -- DON'T compute at poles, blows up dividing by cos(90); UNITS are in (DEG K/sec) C ALSO, ONLY COMPUTE ABOVE 40 km, NO TROPOSPHERIC DDH.... and only include contribution from C small scale gravity waves, i.e., the Kzz term here (Jan. 1997) do 151 ik=1,59 do 151 ij=1,37 dtt1(ij,ik)= rho59(ik)*xkzzt(ij,ik)*((pres59(1)/pres59(ik))**.286) c *(.286/hh*temp1(ij,ik) + dtdz(ij,ik)) 151 continue call derv4(1, dtt1, dt1z, delz, L$1, Z$59, 2) ccc do 152 ik=1,59 ccc do 152 ij=1,37 ccc 152 dt1y(ij,ik) = cosl5(ij)*kyy(ij,ik)*dtdy(ij,ik) ccc ccc call derv4(1, dt1y, dt1yy, dely, L$1, Z$59, 1) do 155 ik=21,59 do 155 ij=2,36 ddh(ij,ik) = -((pres59(ik)/pres59(1))**.286)/rho59(ik)*dt1z(ij,ik) ccc ddh(ij,ik) = -dt1yy(ij,ik)/cosl5(ij) ! old version includes DDH from Kyy, Kzz ccc c - ((pres59(ik)/pres59(1))**.286)/rho59(ik)*dt1z(ij,ik) 155 continue CALL SMOOTH5(ddh, ddhs, 37, 59) CALL SMOOTH5(heat1, heat1s, 37, 59) c smooth ddh and heat1 (Joan's heat rates) terms 1X each C C TOTAL HEATING ;; qtot, heat1s, ddhs, lh1 are in K/sec C do 175 ik=1,59 do 175 ij=1,37 175 qtot(ij,ik) = heat1s(ij,ik) - ddhs(ij,ik) + lh1(ij,ik) c smooth total heating rates 2X cc CALL SMOOTH5(qtot, qtots, 37, 59) cc CALL SMOOTH5(qtots, qtots2, 37, 59) c C C TOTAL MOMENTUM FORCING C do 180 ik=1,59 do 180 ij=1,37 C kru(ij,ik) = kry(ij,ik)*ubar1(ij,ik) kru(ij,ik) = 0.0 ! set Rayleigh Friction = 0 180 fxall(ij,ik) = -kru(ij,ik) + fxtt(ij,ik) + qv(ij,ik) + fk(ij,ik) C Compute derivatives of Forcing C call derv4(1, qtot, dhdy, dely, L$1, Z$59, 1) call derv4(1, fxall, dfxdz, delz, L$1, Z$59, 2) C C ff1 = ff0 + 2*Ubar*tan(xlat5)/a ; ff2 = ff0 + ubar*tan(xlat5)/a - du/dy ;; ff0 = f C df1dy is in analytic form after expanding in spherical coords, C do 190 ik=1,59 do 190 ij=2,36 ff1(ij,ik) = ff0(ij) + 2.*ubar1(ij,ik)*tana(ij) ff2(ij,ik) = ff0(ij) + ubar1(ij,ik)*tana(ij) - dudy(ij,ik) df1dy(ij,ik) = 2.*om1*cosl5(ij)/rad + 2.*tana(ij)*dudy(ij,ik) c + 2.*ubar1(ij,ik)/(rad*rad*cosl5(ij)*cosl5(ij)) 190 continue C set f's, dfdy at poles so they don't blow up, doesn't matter since they're not used in X* calculation do 191 ik=1,59 ff1(1,ik) = ff0(1) ff1(37,ik) = ff0(37) ff2(1,ik) = ff0(1) ff2(37,ik) = ff0(37) df1dy(1,ik) = 0.0 191 df1dy(37,ik) = 0.0 C C NOW SET UP COEFFICIENTS FOR STREAMFUNCTION EQUATION, note: forcing term includes cos(xlat5) factor C do 200 ik=1,59 do 200 ij=1,37 czz(ij,ik) = ff1(ij,ik)*ff2(ij,ik) cz(ij,ik) = -ff1(ij,ik)*ff2(ij,ik)/hh + df1dy(ij,ik)*dudz(ij,ik) czy(ij,ik) = 2.*ff1(ij,ik)*dudz(ij,ik) cy(ij,ik) = (xn2(ij,ik) + rr1/hh*dtdz(ij,ik))*tana(ij) - c ff1(ij,ik)/hh*dudz(ij,ik) cyy(ij,ik) = xn2(ij,ik) + rr1/hh*dtdz(ij,ik) c0(ij,ik) = -df1dy(ij,ik)*dudz(ij,ik)/hh cff(ij,ik) = (ff1(ij,ik)*dfxdz(ij,ik) + rr1/hh*dhdy(ij,ik))* c cosl5(ij) 200 continue c reset cy at poles (prop. to tan(lat)), but its not used below anyway do 203 ik=1,59 cy(1,ik) = cy(2,ik) 203 cy(37,ik) = cy(36,ik) c smooth forcing term since its based on derivatives CALL SMOOTH5(cff, cffs, 37, 59) c damp out forcing for levels 55-58 to be 0 at top (level 59) do 801 ij=1,37 ! doesn't help smooth to set forcing to zero at levs 59 AND 58 cffs(ij,59) = 0.0 do 801 ik=55,58 801 cffs(ij,ik) = -cffs(ij,54)/5.*(ik-59) c C C BOUNDARY CONDITIONS, RESET FOR ALL TIME STEPS: Set X* = 0 at poles, TOP and bottom for rigid walls C use G92 eq(9) for LEVEL 2 FIRST GUESS, interpolate between +-15 deg, C USE ONLY BBC = 0 AT LEVEL 1, LEVEL 2 IS specified c C DO 220 ij=2,36 C sum=0.0 C do 221 ik=2,59 C221 sum = sum + rho59(ik)* C c (-kru(ij,ik) + fxtt(ij,ik) + qv(ij,ik) + fk(ij,ik))*delz C C if (ff0(ij) .ne. 0.0) bbc2(ij) = -cosl5(ij)/(ff0(ij)*rho59(2))*sum C220 CONTINUE C C do 222 ij=17,21 C 222 bbc2(ij) = (bbc2(22) - bbc2(16))/6.*(ij-16.) + bbc2(16) C c do 223 ij=1,37 c 223 bbc2(ij) = 0.d0 C Initialize chi field, start with zero's everywhere except bottom boundary do 225 itt=1,2 do 225 ik=1,59 do 225 ij=1,37 225 chi0(ij,ik,itt) = 0.d0 C ITERATION PROCEDURE, iterate until change from between iterations is <1% at all grid points C Set boundary conditionS for first step, bottom, top = 0.0, specified for lev 2 DO 1000 nt1 = 1,iterc do 230 ij=1,37 chi0(ij,59,1) = 0.d0 chi0(ij,1,1) = 0.d0 230 chi0(ij,2,1) = DBLE(bbcadj(ij,im10)) c DON'T solve for X* at poles, top, or bottom, or level 2, these are pre-set do 400 ik=3,58 do 400 ij=2,36 sf1=chi0(ij-1,ik-1,2)*(czy(ij,ik)*dely*delz) sf2=chi0(ij,ik-1,2)*(4.*czz(ij,ik)*dely2 -2.*cz(ij,ik)*dely2*delz) sf3=chi0(ij+1,ik-1,2)*(-czy(ij,ik)*dely*delz) sf4=chi0(ij-1,ik,2)*(4.*cyy(ij,ik)*delz2 -2.*cy(ij,ik)*dely*delz2) sf6=chi0(ij+1,ik,1)*(2.*dely*delz2*cy(ij,ik) +4.*delz2*cyy(ij,ik)) sf7=chi0(ij-1,ik+1,1)*(-dely*delz*czy(ij,ik)) sf8=chi0(ij,ik+1,1)*(4.*dely2*czz(ij,ik) +2.*dely2*delz*cz(ij,ik)) sf9=chi0(ij+1,ik+1,1)*(dely*delz*czy(ij,ik)) sf0 = 4.*dely2*delz2*cffs(ij,ik) sf5 = 8.*dely2*czz(ij,ik) + 8.*delz2*cyy(ij,ik) c - 4.*dely2*delz2*c0(ij,ik) chi0(ij,ik,2) = chi0(ij,ik,1) + srr/sf5*(sf1 + sf2 + sf3 + sf4 c + sf6 + sf7 + sf8 + sf9 - sf0 - sf5*chi0(ij,ik,1)) 400 continue C C SET Boundaries = 0 for all time steps do 260 ik=1,59 chi0(1,ik,2) = 0.0 260 chi0(37,ik,2) = 0.0 do 270 ij=1,37 cc chi0(ij,59,2) = chi0(ij,58,2) chi0(ij,59,2) = 0.d0 chi0(ij,1,2) = 0.d0 270 chi0(ij,2,2) = DBLE(bbcadj(ij,im10)) C Now check difference between iterations, after 100 iterations, first initialize each time C if maximum difference is <0.1%, end iteration loop, and goto w, v calculation if (nt1 .gt. 100) then do 280 ik=2,58 do 280 ij=2,36 diffc(ij,ik) = 0.0d0 if (chi0(ij,ik,1) .ne. 0.0) c diffc(ij,ik)= 100.*(chi0(ij,ik,2)-chi0(ij,ik,1))/chi0(ij,ik,1) 280 continue dtest = 0.0 do 285 ik=2,58 do 285 ij=2,36 adiffc = DABS(diffc(ij,ik)) 285 dtest = DMAX1(adiffc, dtest) if (dtest .lt. 0.1) go to 99 endif c if continuing, update chi0 for next iteration do 295 ik=1,59 do 295 ij=1,37 295 chi0(ij,ik,1) = chi0(ij,ik,2) 1000 CONTINUE 99 do 305 ik=1,59 do 305 ij=1,37 305 chif(ij,ik) = chi0(ij,ik,2) c Compute w*, v* from stream function (A-grid), just for diagnostics walla, valla(37,59,36) call derv4(1, chif, dcdz, delz, L$1, Z$59, 2) do 307 ik=1,59 do 307 ij=2,36 valla(ij,ik,im10) = c -(dcdz(ij,ik) - chif(ij,ik)/hh)/cosl5(ij)*100. ! also convert to cm/sec 307 continue do 309 ik=1,59 valla(1,ik,im10) = 0.0 309 valla(37,ik,im10) = 0.0 c Now Compute w* from stream function, A-grid = walla(37,59,36)=real*4 , wchi=(37,59)=real*8 c call derv4(1, chif, dcdy, dely, L$1, Z$59, 1) do 310 ik=2,59 do 310 ij=2,36 310 wchi(ij,ik) = dcdy(ij,ik)/cosl5(ij)*100. ! convert to cm/sec here do 311 ij=1,37 wchi(ij,1) = 0.0d0 ! top and bottom w* = 0, 311 wchi(ij,59) = 0.0d0 do 312 ik=1,59 wchi(1,ik) = wchi(2,ik) ! just set w*=85deg at poles, just for diagnostics 312 wchi(37,ik) = wchi(36,ik) do 313 ik=1,59 do 313 ij=1,37 313 walla(ij,ik,im10) = wchi(ij,ik) c ! Now w* on C-grid, Reduce to 18 latitudes do 315 ik=1,59 do 315 ij=1,18 315 w1(ij,ik) = wchi(ij+ij,ik) c c Ramp down w* at the top 5 levels, set w*=0 at level 59 (Top) bottom is already zero from above c c DO 250 ij=1,18 c w1(ij,55) = .8*w1(ij,55) c w1(ij,56) = .6*w1(ij,56) c w1(ij,57) = .4*w1(ij,57) c w1(ij,58) = .2*w1(ij,58) c250 w1(ij,59) = 0.0d0 c c subroutine vwmass does mass adjustment at each p-level, and computes v* by simple mass continuity c NOTE: THIS IS DONE HERE JUST FOR THE OUTPUT FILES (WALL, VALL). VWMASS is then called again C every day in SETDAILY, after w* is interpolated to daily values C CALL VWMASS print *, 'im10, streamf iterations = ', im10, nt1 cc if (mod(im10-2,9) .eq. 0.0) then cc do 199 ik=1,59 cc write(6,299) ik, WBAR(1,ik), WBAR(2,ik), w1(13,ik), w(13,ik) cc 299 format(1x, i5, 2x, 1P4D14.4) cc 199 continue cc print *, ' ' cc endif c now load/reduce arrays for use in model, need to do WALL and VALL separately, w(18,59), v(19,58) c VALL here is just used for output, it's updated every day in model do 600 ik=1,59 do 600 ij=1,18 600 wall(ij,ik,im10) = w(ij,ik) do 610 ik=1,58 do 610 ij=1,19 610 vall(ij,ik,im10) = v(ij,ik) c convert Kyy, Kzz to cm2/sec do 501 ik=1,59 do 501 ij=1,18 kzzall(ij,ik,im10) = xkzzt(ij+ij,ik)*1.d4 501 tempall(ij,ik,im10) = temp10(ij+ij,ik,im10) c For reducing Kyy's to 18 latitudes: Since we're losing resolution here, C assume 35N,S = max value over 30-40N,S interval, to get subtropical surf zone, then c 45N,S = 45-50N,S avg, 55N,S=55/60 avg, 65N,S=65/70 avg, 75N,S=75/80 avg, 85N,S and 25S-25N are as is c This is better overall than just skipping every other latitude, kyy = (37,59) c do 510 ik=1,59 do 505 ij=7,12 505 kyyall(ij,ik,im10) = kyy(ij+ij,ik)*1.d4 kyyall(1,ik,im10) = kyy(2,ik)*1.d4 kyyall(18,ik,im10) = kyy(36,ik)*1.d4 do 506 ij=2,5 506 kyyall(ij,ik,im10) = (kyy(ij+ij-1,ik) + kyy(ij+ij,ik))/2.*1.d4 do 504 ij=14,17 504 kyyall(ij,ik,im10) = (kyy(ij+ij,ik) + kyy(ij+ij+1,ik))/2.*1.d4 kk1 = DMAX1(kyy(11,ik), kyy(12,ik), kyy(13,ik)) kyyall(6,ik,im10) = kk1*1.d4 kk2 = DMAX1(kyy(25,ik), kyy(26,ik), kyy(27,ik)) kyyall(13,ik,im10) = kk2*1.d4 510 CONTINUE c now load other dynamics arrays NOT USED IN MODEL, into full annual arrays, don't reduce to 18 latitudes do 507 ik=1,59 do 507 ij=1,37 ubarall(ij,ik,im10) = ubar1(ij,ik) kzzall37(ij,ik,im10) = xkzzt(ij,ik) fxtall(ij,ik,im10) = fxtt(ij,ik) fkall(ij,ik,im10) = fk(ij,ik) qvall(ij,ik,im10) = qv(ij,ik) kruall(ij,ik,im10) = kru(ij,ik) heatall(ij,ik,im10) = heat1s(ij,ik) ddhall(ij,ik,im10) = ddhs(ij,ik) lhall(ij,ik,im10) = lh1(ij,ik) qtotall(ij,ik,im10) = qtot(ij,ik) cffall(ij,ik,im10) = cffs(ij,ik) chiall(ij,ik,im10) = chif(ij,ik) 507 continue SAVE RETURN END cdbc The following subroutine is copied from the 2d model written by cdbc Mark Schoeberl. (copied 10/22/91) It has been modified so cdbc that boundaries are not set to zero. -- TAKEN FROM OLD KYY ROUTINE SUBROUTINE DERV4(TYP,XIN,OUT,DEL,N,M,SW) C DO A FOURTH OR SECOND ORDER FIRST DERIVATIVE C C IF TYP = 1 THEN SECOND ORDER C IF TYP = 2 THEN FOURTH ORDER C C IF SW=1 THEN INNER INDEX C IF SW=2 THEN OUTER INDEX REAL*8 xin(N,M), out(N,M), del REAL*8 ay1, ay2, ay3 INTEGER N,M INTEGER TYP,SW GO TO (1,2),TYP STOP ' WRONG TYP IN DERV4' 2 AY1=1./(12.*DEL) AY2=2./(3.*DEL) AY3=1./(2.*DEL) GOTO 8 1 AY1=0. AY2=1./(2.*DEL) AY3=AY2 8 GO TO (10,20),SW STOP 'ILLEGAL CALL TO DERV4' 10 N1=N-1 N2=N-2 C - F. Vitt correction put in 5/19/94 DO 51 IK=1,M DO 50 IJ=3,N2 OUT(IJ,IK)=(-AY1*XIN(IJ+2,IK)+AY2*XIN(IJ+1,IK)-AY2*XIN(IJ-1,IK) 1 +AY1*XIN(IJ-2,IK)) 50 CONTINUE OUT(2,IK)=(XIN(3,IK)-XIN(1,IK))*AY3 OUT(n1,IK)=(XIN(n,IK)-XIN(n2,IK))*AY3 OUT(1,IK)=(xin(2,ik)-xin(1,ik))/del OUT(N,IK)=(xin(n,ik)-xin(n1,ik))/del 51 CONTINUE RETURN 20 M1=M-1 M2=M-2 DO 61 IJ=1,N DO 60 IK=3,M2 OUT(IJ,IK)=(-AY1*XIN(IJ,IK+2)+AY2*XIN(IJ,IK+1)-AY2*XIN(IJ,IK-1) 1+AY1*XIN(IJ,IK-2)) 60 CONTINUE OUT(IJ,2)=(XIN(IJ,3)-XIN(IJ,1))*AY3 OUT(IJ,m1)=(XIN(IJ,m)-XIN(IJ,m2))*AY3 OUT(IJ,1)=(xin(ij,2)-xin(ij,1))/del OUT(IJ,M)=(xin(ij,m)-xin(ij,m1))/del 61 CONTINUE RETURN END