subroutine setvbc(j) c c======================================================================= c === c Compute vertical boundary conditions for momentum (CLINIC) and === c tracer (TRACER) equations. === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #include #include #ifdef forcing # include #endif #if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse # include #endif #include #include #include #if defined resetjulian # include #endif #if defined ext_tide & defined mixtide # include #endif c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c logical first integer i,j,m,kz #ifdef forcing integer it1,it2 FLOAT * dstarf,dt1,dt2,dt3,fx,fy,snfscl,sqfscl,swfscl #endif FLOAT * uvmag c save first #ifdef forcing save dstarf,snfscl,swfscl,sqfscl #endif #ifdef bioDuse FLOAT * wsnkphy2 #endif data first /.true./ c c======================================================================= c Begin executable code. c======================================================================= #ifdef forcing c c----------------------------------------------------------------------- c On first call, open input forcing NetCDF and read in first set of c forcing fields. Compute factors to convert fluxes. c----------------------------------------------------------------------- c if(first) then first=.false. call readvbc sqfscl=c1e6/(rho0*cp) swfscl=sec2day snfscl=sec2day/c10 #ifndef resetjulian dstarf=dstart #else dstarf=d0start #endif endif c c----------------------------------------------------------------------- c If appropriate, read in next set of forcing fields. c----------------------------------------------------------------------- c call readvbc #endif c c----------------------------------------------------------------------- c Set surface momentum flux: wind stress (dynes/cm^2). c----------------------------------------------------------------------- #ifdef forcing c c Time-interpolate surface momentum flux from gridded or point data. c if(lsmflux) then it1=1-itsmf it2=itsmf dt1=(tsmf(it2)-dstarf)*day2sec-ttsec dt2=ttsec-(tsmf(it1)-dstarf)*day2sec dt3=dt1+dt2 if(lsmfgrd) then do 10 i=1,imt smf(i,1)=(smflux(1,i,j,it1)*dt1+smflux(1,i,j,it2)*dt2)/dt3 smf(i,2)=(smflux(2,i,j,it1)*dt1+smflux(2,i,j,it2)*dt2)/dt3 wsx(i)=smf(i,1) wsy(i)=smf(i,2) 10 continue else fx=(smfx(1,it1)*dt1+smfx(1,it2)*dt2)/dt3 fy=(smfx(2,it1)*dt1+smfx(2,it2)*dt2)/dt3 do 20 i=1,imt smf(i,1)=fx smf(i,2)=fy wsx(i)=smf(i,1) wsy(i)=smf(i,2) 20 continue endif else do 30 i=1,imt smf(i,1)=c0 smf(i,2)=c0 30 continue endif #else c c Set no-flux condition of momentum at the surface. c do 30 i=1,imt smf(i,1)=c0 smf(i,2)=c0 30 continue #endif c c----------------------------------------------------------------------- c Set bottom momentum flux: bottom stress. c----------------------------------------------------------------------- c do 40 i=1,imt kz=kmu(i) if(kz.ne.0) then #ifndef bottom # if defined ext_tide & defined mixtide uvmag=sqrt(ub(i,kz)**2+vb(i,kz)**2+p5*btide(i)) # else uvmag=sqrt(ub(i,kz)**2+vb(i,kz)**2) # endif bmf(i,1)=cdbot*ub(i,kz)*uvmag bmf(i,2)=cdbot*vb(i,kz)*uvmag #else # if defined ext_tide & defined mixtide uvmag=sqrt(u(i,kz)**2+v(i,kz)**2+p5*btide(i)) # else uvmag=sqrt(u(i,kz)**2+v(i,kz)**2) # endif cmf(i)=cdbot*uvmag bmf(i,1)=cmf(i)*ub(i,kz) bmf(i,2)=cmf(i)*vb(i,kz) #endif else #ifdef bottom cmf(i)=c0 #endif bmf(i,1)=c0 bmf(i,2)=c0 endif 40 continue c c----------------------------------------------------------------------- c Set surface (net) heat flux (convert from W/m2 to degC cm/sec). c 1 degC cm sec-1 ~= 1 cal cm-2 sec-1 c----------------------------------------------------------------------- #ifdef forcing c c Time-interpolate surface heat flux from gridded or point data. c if(lshflux) then it1=1-itshf it2=itshf dt1=(tshf(it2)-dstarf)*day2sec-ttsec dt2=ttsec-(tshf(it1)-dstarf)*day2sec dt3=dt1+dt2 if(lshfgrd) then do 50 i=1,imt stf(i,1)=sqfscl*(shflux(i,j,it1)*dt1+shflux(i,j,it2)*dt2)/ * dt3 stfp(i,1)=sqfscl*(shflux(i,j+1,it1)*dt1+shflux(i,j+1,it2)* * dt2)/dt3 50 continue else fx=sqfscl*(shfx(it1)*dt1+shfx(it2)*dt2)/dt3 do 60 i=1,imt stf(i,1)=fx stfp(i,1)=fx 60 continue endif else do 70 i=1,imt stf(i,1)=c0 stfp(i,1)=c0 70 continue endif #else c c Set no-flux condition of heat at the surface. c do 70 i=1,imt stf(i,1)=c0 stfp(i,1)=c0 70 continue #endif c c----------------------------------------------------------------------- c Set surface water flux: E-P (convert from cm/day to ppt cm/s). c----------------------------------------------------------------------- #ifdef forcing c c Time-interpolate E-P from gridded or point data. c if(lswflux) then it1=1-itswf it2=itswf dt1=(tswf(it2)-dstarf)*day2sec-ttsec dt2=ttsec-(tswf(it1)-dstarf)*day2sec dt3=dt1+dt2 if(lswfgrd) then do 80 i=1,imt fx=swfscl*(swflux(i,j,it1)*dt1+swflux(i,j,it2)*dt2)/dt3 stf(i,2)=(tb(i,1,2)+smean)*fx fx=swfscl*(swflux(i,j+1,it1)*dt1+swflux(i,j+1,it2)*dt2)/dt3 stfp(i,2)=(tbp(i,1,2)+smean)*fx 80 continue else fx=swfscl*(swfx(it1)*dt1+swfx(it2)*dt2)/dt3 do 90 i=1,imt stf(i,2)=(tb(i,1,2)+smean)*fx stfp(i,2)=(tbp(i,1,2)+smean)*fx 90 continue endif else do 100 i=1,imt stf(i,2)=c0 stfp(i,2)=c0 100 continue endif #else c c Set no-flux condition of salt at the surface. c do 100 i=1,imt stf(i,2)=c0 stfp(i,2)=c0 100 continue #endif #if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse c c----------------------------------------------------------------------- c Set surface fluxes for biological tracers. c----------------------------------------------------------------------- c # ifdef forcing c c Set dilution flux due to surface water flux, if any. c do 110 m=3,nt do 110 i=1,imt fx=stf(i,2)/(tb(i,1,2)+smean) stf(i,m)=tb(i,1,m)*fx 110 continue c c Time-interpolate nitrogen flux from gridded or point data. c if(lsnflux) then it1=1-itsnf it2=itsnf dt1=(tsnf(it2)-dstarf)*day2sec-ttsec dt2=ttsec-(tsnf(it1)-dstarf)*day2sec dt3=dt1+dt2 if(lsnfgrd) then do 120 i=1,imt fx=snfscl*(snflux(i,j,it1)*dt1+snflux(i,j,it2)*dt2)/dt3 stf(i,ino3)=stf(i,ino3)+fx 120 continue else fx=snfscl*(snfx(it1)*dt1+snfx(it2)*dt2)/dt3 do 130 i=1,imt stf(i,ino3)=stf(i,ino3)+fx 130 continue endif endif # else c c Set no-flux condition of biological tracers at the surface. c do 140 m=3,nt do 140 i=1,imt stf(i,m)=c0 140 continue # endif #endif c c----------------------------------------------------------------------- c Set bottom tracer flux. c----------------------------------------------------------------------- c c Set no-flux condition of tracers at the bottom. c do 150 m=1,nt do 150 i=1,imt btf(i,m)=c0 #ifdef bioAnder if(m.eq.iphy) btf(i,m)=wsnkphy*tb(i,km,iphy) if(m.eq.idet) btf(i,m)=wsnkdet*tb(i,km,idet) if(m.eq.ino3) btf(i,m)=-fracrmn*(wsnkphy*tb(i,km,iphy)+ * wsnkdet*tb(i,km,idet)) #endif #ifdef bioDuse cjad why not use kz=kmt(i) here instead of km?? (i.e. tracer.F) if(wsnkphy.lt.c0) then if(tb(i,km,ichl).gt.c0) then wsnkphy2=c2*(chl2rad*tb(i,km,ichl)**chl2rade)* * grav*cm2m* * (rhophy-rhos(i,km) # ifdef rmdenbar * -rhobar(i,j,km) # endif * )/(c9*visch2o) else wsnkphy2=c0 endif else wsnkphy2=wsnkphy endif if(biopos.eq.2) then if(m.eq.ichl) btf(i,m)=wsnkphy2*max(c0,tb(i,km,ichl)) if(m.eq.iqn3) btf(i,m)=wsnkphy2*max(c0,tb(i,km,iqn3)) if(m.eq.iqn4) btf(i,m)=wsnkphy2*max(c0,tb(i,km,iqn4)) if(m.eq.idet) btf(i,m)=wsnkdet*max(c0,tb(i,km,idet)) if(m.eq.ino3) btf(i,m)=-fracrmn*(wsnkphy2* * (max(c0,tb(i,km,iqn3)) * +max(c0,tb(i,km,iqn4))) * +wsnkdet*max(c0,tb(i,km,idet))) else if(m.eq.ichl) btf(i,m)=wsnkphy2*tb(i,km,ichl) if(m.eq.iqn3) btf(i,m)=wsnkphy2*tb(i,km,iqn3) if(m.eq.iqn4) btf(i,m)=wsnkphy2*tb(i,km,iqn4) if(m.eq.idet) btf(i,m)=wsnkdet*tb(i,km,idet) if(m.eq.ino3) btf(i,m)=-fracrmn*(wsnkphy2*(tb(i,km,iqn3) * +tb(i,km,iqn4)) * +wsnkdet*tb(i,km,idet)) endif #endif #ifdef bioFasham if(m.eq.ipon) btf(i,m)=wsnkpon*tb(i,km,ipon) if(m.eq.inh4) btf(i,m)=-wsnkpon*tb(i,km,ipon) #endif 150 continue #ifdef forcing c c----------------------------------------------------------------------- c Set shortwave radiation (convert from W/m2 to degC cm/sec). c----------------------------------------------------------------------- c c Time-interpolate shortwave radiation from gridded or point data. c if(lswrad) then it1=1-itswr it2=itswr dt1=(tswr(it2)-dstarf)*day2sec-ttsec dt2=ttsec-(tswr(it1)-dstarf)*day2sec dt3=dt1+dt2 if(lswrgrd) then do 160 i=1,imt srf(i)=sqfscl*(swrad(i,j,it1)*dt1+swrad(i,j,it2)*dt2)/dt3 srfp(i)=sqfscl*(swrad(i,j+1,it1)*dt1+swrad(i,j+1,it2)*dt2)/ * dt3 160 continue else fx=sqfscl*(swr(it1)*dt1+swr(it2)*dt2)/dt3 do 170 i=1,imt srf(i)=fx srfp(i)=fx 170 continue endif else do 180 i=1,imt srf(i)=c0 srfp(i)=c0 180 continue endif #else c c Set no-flux condition at the surface. c do 190 i=1,imt srf(i)=c0 srfp(i)=c0 190 continue #endif return end