c ------------------------------------------------------------------------- subroutine hflx_s89(tenso,npt,iox,t,tclim,cloudy,y,taux,tauy,q,qr,qb,tp) c ------------------------------------------------------------------------- c update the model heat flux using the 1989 Seager et al. formulation. c t = (input) model sst in degrees c. c cloudy = (input) cloud cover as a fraction (0 to 1). c y = (input) model grid latitudes in degrees. c q = (output) heat flux. dimension iox(1),t(1),tclim(1),cloudy(1),y(1),taux(1),tauy(1), * q(1),qr(1), qb(npt,1),tp(nyp,1) complex CPI18,EXLAT,EXPHI,EX2PHI common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common /new_hfxevp/ QCON, rlx_time, solr_gamma, TATM, SATM * ,trans_coef_sst,trans_coef_sss save jmon parameter (RTODEG = 180./3.14159265) data TAUA/.017/,TAUCON/10300./,UVMIN/4.0/ data cpi18 /(0.98481,-0.17365)/ c......for the tropic pacific: data al,acld /0.06,0.75/ data aalpha,rh,asst,ct0 /0.0020,0.78,1.667,-5.0/ c......for the tropical atlantic: c data al,acld /0.06,0.66/ c data aalpha,rh,asst,ct0 /0.0019,0.735,1.8,-5.0/ c.....parameter definitions: c rjuldy = julian day (0 = 1 jan). c uvmin = minimum wind speed in m/s. c al = albedo. c acld = cloud coefficient. c aalpha = alpha coefficient. c rh = pseudo relative humidity. c asst = sst coefficient. c cto = sst coefficient times offset. c cpi18 = cexp( (0,1)*(-pi/18.0) ). c c.....set constants used to find heat flux. a1mrh = 1.0 - rh a1mal = 1.0 - al c.....get the current model time in julian days: call DayOfYear(tenso, juld, july) call enso2date(tenso, id,imon,iy) rjuldy = juld-1 pha = 6.28318/real(july) phi = pha*(rjuldy - 21.0) exphi = cexp((0,1)*phi) cosphi = real(exphi) sinphi = aimag(exphi) ex2phi = exphi**2 cos2phi = real(ex2phi) sin2phi = aimag(ex2phi) arg = 23.45*sin((rjuldy-82.0)*pha) c the formulae from Weare 1980 are converted to w/m**2 and radians. c e.g. a1 from Weare is now a1=9.63 + 192.44*cos(rlat+90) and c in the calculatiion for the noon solar angle c alpha = arcsin(cos(rlat)*cos(arg1) + sin(rlat)*sin(arg1)) c = arcsin(cos(rlat - arg1)) c = arcsin(sin(rlat - arg1 + 90.)). c c.....compute the latitude dependent coefficients. do j = 1, nyp rlat = y(j)/rtodeg exlat = cexp((0,1)*rlat) cosrlat = real(exlat) sinrlat = aimag(exlat) c........set the coefficients for the clear sky solar radiation function. a0 = -15.82 + 326.87*cosrlat a1 = 9.63 - 192.44*sinrlat b1 = -3.27 + 108.7*sinrlat a2 = -0.64 -7.8*real(exlat**2) b2 = -0.5 + 14.42*real((exlat**2)*cpi18) c........tp(j,1) = q0 ! find q0, the clear sky radition. tp(j,1) = a0 + a1*cosphi + b1*sinphi + a2*cos2phi + b2*sin2phi c........find alpha, the noon solar altitude or angle. alpha = arg - y(j) + 90. tp(j,2) = alpha enddo c c this if statement is equivalent to arcsin(sin(x)) for the range c of alpha encountered. do 20 j = 1, nyp 20 if(tp(j,2) .gt. 90.0) tp(j,2) = 180. - tp(j,2) taud2 = sqrt(TAUCON/TAUA) qcon_inv = 1. / qcon qcon_gam = SOLR_GAMMA / qcon cnst1 = 2500000.0 / 461.0 cnst2 = 1.0 / 273.15 cnst3 = 0.622*6.11 cnst4 = 4.59373*a1mrh c.....compute the heat flux. do k = 1, npt j = (iox(k)-1)/nxp + 1 i = iox(k) - (j-1)*nxp c........set a minimum wind speed of UVMIN(m/s) to avoid underestimate c of wind magnitude due to use of monthly averages. tc = amax1(taud2*(taux(k)**2 + tauy(k)**2)**.25, UVMIN) exparg = cnst1*(cnst2 - 1.0/(t(k)+273.15)) qs = cnst3 * exp(exparg) c = a1mal*(1.0 - acld*cloudy(k) + aalpha*tp(j,2)) qsolr = c*tp(j,1) qb(k,1) = qsolr qr(k) = qcon_gam * qsolr qb(k,2) = cnst4 * tc * qs qb(k,3) = asst*t(k) + ct0 qb(k,5) = 30.*(tclim(k)-t(k)) qtot = qsolr - qb(k,2) - qb(k,3) q(k) = qcon_inv * qtot - qr(k) c if ( imon.ne.jmon .and. j.eq.nyp/2 .and. i.eq.nxp/2) c * write (2, '(1hq,i2,5(1pg12.4))') imon,qtot,qb(k,1), c * qb(k,2), qb(k,3), qb(k,5) enddo jmon = imon return end subroutine hflx_s94 (npt, t, taux,tauy,sst,cld,solr, q,qr,qb) c------------------------------------------------------------------------- c according to R.Seager & B.Blumental, December 1994, Journal of Climate. c q = (output) heat flux = f(t,cld,solr) - gamma*solr dimension t(1),taux(1),tauy(1),sst(1),cld(1),solr(1),q(1),qr(1),qb(npt,1) common /new_hfxevp/ QCON, rlx_time, solr_gamma, TATM, SATM * ,trans_coef_sst,trans_coef_sss data TAUA/.017/,TAUCON/10300./,UVMIN/4.0/ data AL,ALATENT,CE,RHOA/0.06,2.5e06,1.24e-03,1.225/ data EPS,SIG,TSTA,D /0.97,5.6696e-08,0.71,0.78/ data CP /1004/ cnst1 = 0.622*6.112/1000. cnst2 = RHOA*CE*ALATENT*(1.- D) cnst3 = RHOA*CE*CP*TSTA cnst4 = D*1000./0.622 cnst5 = EPS*SIG cnst6 = 4.*tsta taud2 = sqrt(TAUCON/TAUA) qcon_inv = 1./QCON qcon_gam = SOLR_GAMMA / QCON do i = 1, npt c........Solar Radiation (for output only): qsol = solr(i) qb(i,1) = qsol c........Latent heat: wnsp = amax1(taud2*(taux(i)**2 + tauy(i)**2)**.25, UVMIN) qsat = cnst1*exp(17.67*t(i)/(t(i)+243.5)) qlh = cnst2 * wnsp * qsat qb(i,2) = qlh c........Sensible heat: qsh = cnst3 * wnsp qb(i,3) = qsh c........Long Wave Back Radiation: ts = t(i) + 273.15 sqeth = sqrt(cnst4 * qsat) ts3 = cnst5*ts*ts*ts cnst = 0.8 if( t(i) .gt. 28.) cnst = 0.4 qlw = ts3*(ts*(1.-cnst*cld(i)**2)*(0.417-0.0486*sqeth) + cnst6) qb(i,4) = qlw c........Heat Flux Deficit (for output only): qb(i,5) = 30.*(sst(i)-t(i)) c........The Total Heat Flux at the Surface: qr(i) = qcon_gam * qsol qtot = qsol - qlh - qsh - qlw c c........Since this is in W/m2, we need to convert it to dyn/m3 q(i) = qcon_inv * qtot - qr(i) enddo return end c------------------------------------------------------------------------- subroutine hflx_s94b(npt, t, taux,tauy,sst,cld,solr, q,qr,qb) c------------------------------------------------------------------------- c according to R.Seager & B.Blumental, December 1994, Journal of Climate. c q = (output) heat flux = f(t,cld,solr) - gamma*solr dimension t(1),taux(1),tauy(1),sst(1),cld(1),solr(1),q(1),qr(1),qb(npt,1) common /new_hfxevp/ QCON, rlx_time, solr_gamma, TATM, SATM * ,trans_coef_sst,trans_coef_sss data TAUA/.017/,TAUCON/10300./,UVMIN/4.0/ data AL,ALATENT,CE,RHOA/0.06,2.5e06,1.24e-03,1.225/ data EPS,SIG,TSTA,D /0.97,5.6696e-08,0.71,0.78/ data CP/1004/ cnst1 = 0.622*6.112/1000. cnst2 = RHOA*CE*ALATENT cnst3 = RHOA*CE*CP*TSTA cnst4 = 1000./0.622 cnst5 = EPS*SIG cnst6 = 4.*tsta taud2 = sqrt(TAUCON/TAUA) qcon_inv = 1./QCON qcon_gam = SOLR_GAMMA/QCON do i = 1, npt c........Solar Radiation (for output only): qsol = solr(i) qb(i,1) = qsol c........Latent heat: wnsp = amax1(taud2*(taux(i)**2 + tauy(i)**2)**.25, UVMIN) qsat = cnst1*exp(17.67*t(i)/(t(i)+243.5)) qair = (-9.42 + 0.97*t(i))/1000. qlh = cnst2 * wnsp * (qsat - qair) qb(i,2) = qlh c........Sensible heat: qsh = cnst3 * wnsp qb(i,3) = qsh c........Long Wave Back Radiation: ts = t(i) + 273.15 sqeth = sqrt(cnst4 * qair) ts3 = cnst5*ts*ts*ts cnst = 0.8 if( t(i) .gt. 28.) cnst = 0.4 qlw = ts3*(ts*(1.-cnst*cld(i)**2)*(0.417-0.0486*sqeth) + cnst6) qb(i,4) = qlw c........Heat Flux Deficit (for output only): qb(i,5) = 30.*(sst(i)-t(i)) c........The Total Heat Flux at the Surface: qr(i) = qcon_gam * qsol qtot = qsol - qlh - qsh - qlw c c........Since this is in W/m2, we need to convert it to dyn/m3 q(i) = qcon_inv * qtot - qr(i) enddo return end c--------------------------------------------------------------------------- subroutine init_pbl(npt, NX, NY, xm, ym, iox) c--------------------------------------------------------------------------- implicit real*4(a-h,o-z),integer(i-n) include 'comm_pbl.h' dimension iox(1), lsm1d(1), xm(nx), ym(ny) pointer (p_lsm1d, lsm1d) parameter (TORAD = 3.14159265/180., REARTH = 6378000.) nxy = nx*ny call mem_alloc (p_up, nxy, 2, 'up') call mem_alloc (p_vp, nxy, 2, 'vp') call mem_alloc (p_thv, nxy, 2, 'thv') call mem_alloc (p_the, nxy, 2, 'the') call mem_alloc (p_thve, nxy, 2, 'thve') call mem_alloc (p_thvs, nxy, 2, 'thvs') call mem_alloc (p_pnuxp, nxy, 2, 'pnuxp') call mem_alloc (p_pnuyp, nxy, 2, 'pnuyp') call mem_alloc (p_qe, nxy, 2, 'qe') call mem_alloc (p_qs, nxy, 2, 'qs') call mem_alloc (p_c0, nxy, 2, 'c0') call mem_alloc (p_dx, nxy, 2, 'dx') call mem_alloc (p_dy, ny, 2, 'dy') call mem_alloc (p_lsm, nxy, 1, 'lsm') c determine grid spacing in m ipbl_jsta = 0 ipbl_jend = 0 do j = 1, ny-1 dy(j) = TORAD * REARTH * (ym(j+1)-ym(j)) if (ipbl_jsta .eq. 0 .and. ym(j).gt.pbl_south) ipbl_jsta = j if (ipbl_jend .eq. 0 .and. ym(j).gt.pbl_north) ipbl_jend = j-1 deg2met = TORAD * REARTH * cos(TORAD * 0.5*(ym(j) + ym(j+1))) do i = 1, nx-1 dx(i,j) = deg2met * (xm(i+1) - xm(i)) enddo dx(nx,j) = dx(nx-1,j) enddo j = ny dy(j) = dy(j-1) do i = 1, nx-1 dx(i,j) = deg2met * (xm(i+1) - xm(i)) enddo dx(nx,j) = dx(nx-1,j) if (ym(1) .ge. pbl_south) ipbl_jsta = 2 if (ym(ny) .le. pbl_north) ipbl_jend = ny-1 p_lsm1d = p_lsm do i = 1, npt lsm1d(iox(i)) = i enddo return end c--------------------------------------------------------------------------- subroutine htflux_pbl (npt,nx,ny, iox,xm,ym, sst,cldfr, wspd,u,v, * q,t, rlh,sh,qlw, qa,th) c--------------------------------------------------------------------------- implicit real(a-h,o-z),integer(i-n) include 'comm_pbl.h' dimension xm(nx), ym(ny), iox(1) dimension sst(npt), cldfr(npt), wspd(npt), u(npt), v(npt), * qlw(npt), rlh(npt), sh(npt), * q(nx,ny), t(nx,ny), th(nx,ny), qa(nx,ny) logical FIRST_PBL, ADVEC save FIRST_PBL data FIRST_PBL /.true./ if ( FIRST_PBL ) then print*,'using the Richard-Senya pbl ...' call init_pbl(npt, NX, NY, xm, ym, iox) FIRST_PBL = .false. endif c set model parameters: ADVEC = (ipbl_advec .eq. 1) pnu = pbl_pnu delta = pbl_delta pml = pbl_pml depth = pbl_depth betav = pbl_betav qrad = pbl_grad / 86400. c set constants: r = 287.04 psfc = 100000. rl = 2.5e+6 cp = 1004. rhoa = 1.225 stef = 5.6696e-8 eps = 0.97 jstart = ipbl_jsta jend = ipbl_jend c_new this is just for checking c determine grid spacing in m c slat = ym(1) c conv=2.*3.14/360. c radius=6.37e+6 c dy(1)=radius*(ym(2) - ym(1))*conv c rlat=slat*conv c dx(1,1)=conv*radius*cos(rlat)*(xm(2) - xm(1)) c do 1 i=2,nx c dx(i,1)=conv*radius*cos(rlat)*(xm(i) - xm(i-1)) c 1 continue c do 2 j=2,ny c dyd = ym(j+1) - ym(j) c dy(j)=conv*radius*dyd c rlat=rlat+dyd*conv c dxd = xm(2) - xm(1) c dx(1,j)=conv*radius*cos(rlat)*dxd c do 2 i=2,nx c dxd = xm(i) - xm(i-1) c dx(i,j)=conv*radius*cos(rlat)*dxd c 2 continue c_new end of checking c Two iterations are performed. A smaller exchange coefficient is c used on second iteration if mixed layer is stable. c First find equilibrium values of theta_V and q. These are set to c their observed values over land. do 24 j = 1, ny do 24 i = 1, nx c0(i,j) = .0014 k = lsm(i,j) if (k .ne. 0) then qs(i,j)=.622*6.11*exp(17.67*(1.-243.5/(sst(k)+243.5)))/1000. endif 24 continue iter = 1 itermx = 3 99 continue do 25 j=jstart,jend do 25 i=1,nx if(iter.gt.1 .and. (thv(i,j).gt.thvs(i,j))) then c0(i,j)=.00075 endif k = lsm(i,j) if (k .eq. 0) then c_orig thve(i,j) = (t(i,j)*(psfc/(psfc-.5*pml))**(r/cp))*(1.+.61*q(i,j)) thve(i,j) = t(i,j)*(1.+.61*q(i,j)) qe(i,j) = q(i,j) th(i,j) = t(i,j) else w0 = wspd(k)*pml/depth thvs(i,j) = (sst(k)+273.15) * (1.+.61*qs(i,j)) thve(i,j) = thvs(i,j)+pml*qrad/((1.+betav)*c0(i,j)*w0) qe(i,j) = qs(i,j)/(1.+delta) endif 25 continue c Set equilibrium values to observed at northernmost and southernmost c points. This is required because advection/diffusion cannot be computed c when there is no poleward point. Actual values of air temperature and c air humidity are also set equal to observed values and used in flux c calculation. do 26 i=1,nx do 27 j=1,jstart up(i,j)=0. vp(i,j)=0. pnuxp(i,j)=0. pnuyp(i,j)=0. qe(i,j)=q(i,j) thve(i,j)=(t(i,j)*(1017./(1000.))**(r/cp))*(1.+.61*q(i,j)) qa(i,j)=q(i,j) thv(i,j)=thve(i,j) th(i,j)=t(i,j)*(1017./(1000.))**(r/cp) 27 continue do 26 j=jend,ny up(i,j)=0. vp(i,j)=0. pnuxp(i,j)=0. pnuyp(i,j)=0. qe(i,j)=q(i,j) thve(i,j)=(t(i,j)*(1017./(1000.))**(r/cp))*(1.+.61*q(i,j)) qa(i,j)=q(i,j) thv(i,j)=thve(i,j) th(i,j)=t(i,j)*(1017./(1000.))**(r/cp) 26 continue c Set diffusion and advecting wind speed. Over land both are c set to zero to ensure derived theta_V and q are observed c values. In addition, diffusion is set to zero close to c coastline. do 29 j=jstart,jend do 29 i=1,nx k = lsm(i,j) if ( k .eq. 0 ) then up(i,j) = 0. vp(i,j) = 0. pnuxp(i,j) = 0. pnuyp(i,j) = 0. else w0=wspd(k)*pml/depth ip1=i+1 if(ip1.eq.(nx+1)) ip1=nx ip2=i+2 if(ip2.eq.(nx+2) .or. ip2.eq.(nx+1)) ip2=nx im1=i-1 if(im1.eq.0) im1=1 im2=i-2 if(im2.eq.0 .or. im2.eq.-1) im2=1 jm1=j-1 jm2=j-2 if(jm2.eq.0) jm2=1 jp1=j+1 jp2=j+2 if(jp2.eq.(ny+1)) jp2=ny if( lsm(ip1,j)*lsm(im1,j)*lsm(i,jp1)*lsm(i,jm1)* * lsm(ip2,j)*lsm(im2,j)*lsm(i,jp2)*lsm(i,jm2) .eq. 0 ) then pnuxp(i,j)=0. pnuyp(i,j)=0. else if(i.eq.1 .or. i.eq.nx) then twodx2=dx(i,j)**2. else twodx2=.25*(dx(i,j)+dx(i+1,j))**2. endif pnuxp(i,j)=pnu*pml/((1.+betav)*c0(i,j)*w0*twodx2) pnuyp(i,j)=pnu*pml/((1.+betav)*c0(i,j)*w0*.25* $ (dy(j)+dy(j-1))*(dy(j)+dy(j-1))) endif if ( advec ) then if (u(k) .gt. 0.) then i1=i else i1=i+1 if (i.eq.nx) i1=nx endif up(i,j) = u(k)*pml/((1.+betav)*c0(i,j)*w0*dx(i1,j)) if (v(k) .gt. 0.) then vp(i,j) = v(k)*pml/((1.+betav)*c0(i,j)*w0*dy(j-1)) else vp(i,j) = v(k)*pml/((1.+betav)*c0(i,j)*w0*dy(j)) endif else up(i,j)=0. vp(i,j)=0. endif endif 29 continue c call subroutine that solves for theta_V call adv2Deq1(nx,ny,up,vp,pnuxp,pnuyp,thve,thv) c repeat one time iter=iter+1 if(iter.lt.itermx) goto 99 c set scaled advecting velocities for humidity calculation and c impose no diffusion across continental boundaries do 39 j=jstart,jend do 39 i=1,nx k = lsm(i,j) if (k .ne. 0) then w0=wspd(k)*pml/depth ip1=i+1 if(ip1.eq.(nx+1)) ip1=nx ip2=i+2 if(ip2.eq.(nx+2) .or. ip2.eq.(nx+1)) ip2=nx im1=i-1 if(im1.eq.0) im1=1 im2=i-2 if(im2.eq.0 .or. im2.eq.-1) im2=1 jm1=j-1 jm2=j-2 if(jm2.eq.0) jm2=1 jp1=j+1 jp2=j+2 if(jp2.eq.(ny+1)) jp2=ny if( lsm(ip1,j)*lsm(im1,j)*lsm(i,jp1)*lsm(i,jm1)* * lsm(ip2,j)*lsm(im2,j)*lsm(i,jp2)*lsm(i,jm2) .eq. 0 ) then pnuxp(i,j)=0. pnuyp(i,j)=0. else if(i.eq.1 .or. i.eq.nx) then twodx2=dx(i,j)**2. else twodx2=.25*(dx(i,j)+dx(i+1,j))**2. endif pnuxp(i,j)=pnu*pml/((1.+delta)*c0(i,j)*w0*twodx2) pnuyp(i,j)=pnu*pml/((1.+delta)*c0(i,j)*w0*.25*(dy(j)+dy(j-1))* $ (dy(j)+dy(j-1))) endif if(advec) then if(u(k).gt.0.) then i1=i else i1=i+1 if(i.eq.nx) i1=nx endif up(i,j)=u(k)*pml/((1.+delta)*c0(i,j)*w0*dx(i1,j)) if(v(k).gt.0.) then vp(i,j)=v(k)*pml/((1.+delta)*c0(i,j)*w0*dy(j-1)) else vp(i,j)=v(k)*pml/((1.+delta)*c0(i,j)*w0*dy(j)) endif else up(i,j)=0. vp(i,j)=0. endif endif 39 continue c call solver to derive q call adv2Deq1(nx,ny,up,vp,pnuxp,pnuyp,qe,qa) c calculate theta from theta_V and q c calculate fluxes of sensible and latent heat do j = 1, ny do i = 1, nx k = lsm(i,j) if ( k .ne. 0 ) then th(i,j) = thv(i,j) / (1.+.61*qa(i,j)) co1 = rhoa * c0(i,j) * wspd(k) sstk = sst(k) + 273.15 co2 = sstk - th(i,j) co3 = eps*stef*sstk*sstk*sstk rlh(k) = co1*rl*(qs(i,j)-qa(i,j)) sh(k) = co1*cp*co2 c qlw(k) = eps*stef*(th(i,j)**4.)*(.39-.05*sqrt(qa(i,j)*1000./.622)) c * *(1.-.55*cldfr(k)) + 4.*eps*stef*(th(i,j)**3.)*co2 if ( sstk .gt. 301.15 ) then aout =.4 else aout =.8 endif qlw(k) = co3 * ( sstk * (.417-.0486*sqrt(qa(i,j)*1000./.622)) * * (1. - aout*cldfr(k)*cldfr(k)) + 4.*co2) endif enddo enddo return end c This subroutine computes surface fluxes of latent and sensible heat c in units of W/m^2. The fluxes are computed by a forced advection- c diffusion equation. It solves a equations for the virtual potential c temperature and the air humidity and then inverts the first to get c the air temperature. In both case the balance is one of diffusion, c horizontal advection, surface fluxes and a flux at the mixed layer top. c The mixed layer is a constant depth. c c The model also computes long wave cooling with the Berliand and c Berliand bulk formula (see Seager and Blumenthal, J. Climate, Dec '94 c for example). c c Note added 11/7/94: To date the model has been coupled to an ocean c GCM developed by Ragu Murtugudde, now at GSFC. The results have c been good. Some care is needed at open ocean boundaries it turns out. c In the version as I give it here you will see the computation is done c only for meridional index j=jstart,jend with jstart=25 and jend =ny-1. c This is like putting a boundary in the middle of the southern ocean. c For points poleward of the end points the air humidity and temperature c are set equal to observed values ensuring that values advected in are c realistic. We used ECMWF data at 1000mb. We found that the air-sea c temperature difference given by this data was too large (probably c 'cos the SLP is greater than the lowest ananlysis level of 1000mb) so c we correct it to by a dry adiabaltic lapse rate to an slp of 1017 mb c which corresponded to a reasonable SLP at 40S which is where our ocean c GCM began. Clearly users are free to do whatever they want but c *be cautious*!. c c The limits are to set jstart =2 and jend=ny-1. The end points cannot be c included because of the diffusion operator that would otherwise look c out of array bounds. c c The inputs are: c c sst = array containing the model or observed SST c u = array containing observed low level zonal wind velocity c v = array containing observed low level meridional wind velocity c wspd = array containing observed low level wind speed c lsm = a land sea mask (1=land, 0= ocean) c q = observed low level air humidity (kg/kg) c t = observed low level temperature (K) c cldfr= observed cloud cover c wlat = western latitude of input grid, in degree (e.g. 220.) c slat = southern latitude of input grid, in degrees (e.g. -30.) c dxd = grid spacing in degrees longitude. dxd(i) equals the distance from c the longitude at i-1 to the longitude at i which allows for c uneven grid spacing. c dyd = grid spacing in degrees latitude. dyd(j) equals the distance from c the latitude at j to the latitude at j+1 which allows for c uneven grid spacing. c nx = number of x grid points c ny = number of y grid points c c c The outputs are: c c sh = array containing the sensible heat flux (W/m^2) c rlh = array containing the latent heat flux (W/m^2) c qa = atmospheric mixed layer humidity in kg/kg c th = atmospheric mixed layer potential temperature in K c qlw = longwave radiative heat flux c c Parameters are: c c pnu=diffusivity (m^2/s) c delta - equilibrium q = q0/(1+delta) where q0 is saturation humidity c at the SST c pml=pressure depth (Pa) of the mixed layer c depth=geometric depth of mixed layer = (pml/(rhoa*grav) c qrad=radiative cooling K/s c betav=ratio of downward theta_V flux at mixed layer top to the c surface flux c c0=surface exchange coefficient SUBROUTINE adv2Deq1(nx,ny,UP,VP,NUXP,NUYP,QE,QA) REAL UP(NX,NY),VP(NX,NY), QE(NX,NY) REAL NUXP(NX,NY),NUYP(NX,NY) REAL QA(NX,NY) C variables are dimensioned with X first NXSKP = 1 NYSKP = NX C does X advection C loops over all latitudes IX = 1 DO IY = 1 , NY CALL ADVDIFQ1DX(UP(1,IY),NUXP(1,IY),NX,QE(1,IY),QE(1,IY), * QE(NX,IY),QA(1,IY)) END DO C does Y advection C loops over all longitudes IY = 1 DO IX = 1 , NX C boundary conditions QLEFT = QE(IX,1) QRIGHT = QE(IX,NY) CALL ADVDIFQ1D(VP(IX,1),NUYP(IX,1),NY,QA(IX,1), * QLEFT,QRIGHT,NYSKP,QA(IX,1)) END DO RETURN END SUBROUTINE ADVDIFQ1DX(U2,NU2,NX,QE0,QLEFT,QRIGHT,QA) REAL U2(NX), NU2(NX),QE0(NX), QA(NX), QLEFT, QRIGHT INTEGER NX PARAMETER (MAXDIM=300) REAL*4 AC(MAXDIM),BC(MAXDIM),CC(MAXDIM),QE(MAXDIM) * AUTOMATIC AC,BC,CC,QE IF(NX.GT.MAXDIM)STOP 12 C does inside points NXL1 = NX - 1 DO K = 2 , NXL1 QE(K) = QE0(K) IF(U2(K).GE.0) THEN AC(K) = -U2(K) - NU2(K) BC(K) = 1. + U2(K) + 2*NU2(K) CC(K) = -NU2(K) ELSE C wind blows left AC(K) = -NU2(K) BC(K) = 1. - U2(K) + 2*NU2(K) CC(K) = U2(K) - NU2(K) ENDIF END DO C does left point IF(U2(1).GT.0.0)THEN C boundary condition matters QE(1) = QE0(1) + U2(1)*QLEFT BC(1) = 1. + U2(1) CC(1) = 0.0 ELSE C boundary condition doesn't matter QE(1) = QE0(1) BC(1) = 1. - U2(1) CC(1) = U2(1) ENDIF C does right point IF(U2(NX).GT.0.0)THEN C boundary condition doesn't matter QE(NX) = QE0(NX) BC(NX) = 1. + U2(NX) AC(NX) = -U2(NX) ELSE C boundary condition matters QE(NX) = QE0(NX) - U2(NX)*QRIGHT BC(NX) = 1.0 - U2(NX) AC(NX) = 0.0 ENDIF CALL TRIDAGQA(AC,BC,CC,QE,CC,QE,QA,1,NX) RETURN END C \end{fortran} C \subsection{advq1d} C C Inputs C \begin{variablelist} C \v U U scaled version of $U = (Hu)/(C_e |v| dx)$ C \v NX n_x number of points C \v QE q_e equilibrium (forcing) Q C \v QLEFT {} Qa beyond the left boundary C \v QRIGHT {} Qa beyond the right boundary C \v NXSKP {} integer spacing between consecutive elements of QA C \end{variablelist} C C Output C \begin{variablelist} C \v QA q_a solution for QA C \end{variablelist} C \begin{fortran} SUBROUTINE ADVDIFQ1D(U2,NU2,NX,QE0,QLEFT,QRIGHT,NXSKP,QA) REAL U2(NXSKP,NX), QE0(NXSKP,NX), QA(NXSKP,NX), QLEFT, QRIGHT REAL NU2(NXSKP,NX) INTEGER NX PARAMETER (MAXDIM=300) REAL*4 AC(MAXDIM),BC(MAXDIM),CC(MAXDIM),QE(MAXDIM) * AUTOMATIC AC,BC,CC,QE IF(NX.GT.MAXDIM)STOP 12 C does inside points NXL1 = NX - 1 DO K = 2 , NXL1 QE(K) = QE0(1,K) IF(U2(1,K).GE.0) THEN AC(K) = -U2(1,K) - NU2(1,K) BC(K) = 1. + U2(1,K) + 2*NU2(1,K) CC(K) = - NU2(1,K) ELSE C wind blows left AC(K) = -NU2(1,K) BC(K) = 1. - U2(1,K) + 2*NU2(1,K) CC(K) = U2(1,K)-NU2(1,K) ENDIF END DO C does left point IF(U2(1,1).GT.0.0)THEN C boundary condition matters QE(1) = QE0(1,1) + U2(1,1)*QLEFT BC(1) = 1. + U2(1,1) CC(1) = 0.0 ELSE C boundary condition doesn't matter QE(1) = QE0(1,1) BC(1) = 1. - U2(1,1) CC(1) = U2(1,1) ENDIF C does right point IF(U2(1,NX).GT.0.0)THEN C boundary condition doesn't matter QE(NX) = QE0(1,NX) BC(NX) = 1. + U2(1,NX) AC(NX) = -U2(1,NX) ELSE C boundary condition matters QE(NX) = QE0(1,NX) - U2(1,NX)*QRIGHT BC(NX) = 1.0 - U2(1,NX) AC(NX) = 0.0 ENDIF CALL TRIDAGQA(AC,BC,CC,QE,CC,QE,QA,NXSKP,NX) RETURN END C \end{fortran} C \subsection{tridagqa} C Same as tridag2 in the multiple mode linear equatorial model. Solves C a tridiagonal system. C \begin{fortran} C SUBROUTINE TRIDAGQA(A,B,C,D,E,F,X,NXSKP,N) PARAMETER (JS=1) C C SOLVES THE TRIDIAGONAL SYSTEM C A(I)*X(I-1)+B(I)*X(I)+C(I)*X(I+1)=D(I)...I=JS,N A(JS)=0,C(N)=0 C solves in two passes C pass1 computes E(i), F(i) such that x(i) = E(i)*x(i+1) + F(i) C pass2 computes x(i) from E(i), F(i) C array C can be used as E C array D can be used as F C DIMENSION A(*), B(*), C(*), D(*), E(*), F(*), X(NXSKP,*) C E(JS) = C(JS) / B(JS) F(JS) = D(JS) / B(JS) C JN = JS + 1 DO 10 I = JN, N DN = B(I) - A(I) * E(I-1) E(I) = C(I) / DN F(I) = ( D(I) - A(I) * F(I-1) ) / DN 10 CONTINUE X(1,N) = F(N) XOLD = F(N) IMAX = N-1 DO 20 I = IMAX,JS,-1 XOLD = F(I) - E(I) * XOLD X(1,I) = XOLD 20 CONTINUE RETURN END C \end{fortran} C \end{document}