c ======================================================================= c c AMLICE is consists of subroutines that allow to extent the AML of Seager et al. c to an ice covered ocean. Moreover, the 1D thermodynamic properties of the c ice such as thickness, heat content and concentration c are computed using a Hibler-Oberhuber type formulation c c the main subroutine is called HTFLUXICE c which uses a 1-D thermodynamical ice model called ICETHERMO c and a modified AML called HTFLUXI that calls several solvers c note that one was changed slightly to accomodate arbitray array sizes c c all temperatures are in KELVIN !!!! c c Martin Visbeck and Bob Newton, LDEO, August 22, 1996 c c=========================================================================== subroutine htfluxice(mx,my,nx,ny,lsm,dxd,dyd,tstep,tmonth, + sst,cldfr,wspd,u,v,q,t,rlh,sh,qlw,qsw,pp,qa,th,rh, + sss,qisw,ppi,hice,cice,thice,tsnw,qios,brne,rlhi,shi,qlwi,qswi, + rlc0ice,cpc0ice,qlwice1,qlwice2,iperx) c c====6===1=========2=========3=========4=========5=========6========7==2 c This subroutine computes the surface fluxes of a mixed ice-ocean interface. c The ice model is thermodynamic only, i.e. not ice moves. c The atmospheric boundary layer model is a modifyed version from Seager et al. c and the ice model is constructed after Oberhuber et al. c c All temperatures are in Kelvin. c C Input fields: c c -GRID c mx,my : dimension of arrays c nx,ny : size of the used parts of the array c lsm : land sea mask (1=land, 0=ocean/ice) c dxd,dyd : grid spacing [degree] c slat : southern latitude of grid [degree] c c tstep :delta time [s] c tmonth: month of the year c c -AML+ICE c sst : sea surface temp [K] c cldfr : cloud fraction [0-1] (used for longwave radiation only) c wspd : wind speed [m/s] c u : zonal wind [m/s] c v : meridional wind [m/s] c q : observed humidity [kg/kg] c t : observed temperature [K] c c -ICE c sss : sea surface salinity [psu] c qisw : cloud (but not albedo) corrected incomming short wave radiation [W/m^2] c ppi : precipitation [m/s] C Output fields: c c -ATMOS-OCEAN FLUXES c rlh : latent heat flux [W/m^2] c sh : sensible heat flux [W/m^2] c qlw : long wave net radiation [W/m^2] c qsw : short wave (including albedo) [W/m^2] c pp : precip only over water [m/s] c -ICE-OCEAN FLUXES c qios : ice-ocean flux [W/m^2] c brne : fresh water flux melt/freeze [m/s] c C Additional output for diagnosis c c -AML c th : aml potential temperature [K] c qa : aml humidity [kg/kg] c rh : aml relative humidity [0-1] c -ICE c hice : ice thickness averaged over grid point [m] c cice : ice concentration [0-1] c tsnw : `snow' temperature [K] c thice : ice `heat content' [Km] c -ATMOS-ICE FLUXES c rlhi : latent heat flux ice-atmos [W/m^2] c shi : sensible heat flux ice-atmos [W/m^2] c qlwi : long wave net eadiation ice-atmos [W/m^2] c qswi : short wave (incl. albedo ice-atmos [W/m^2] c c -ATMOS-ICE FLUX parts c rlc0ice : c cpcoice : c qlwice1 : c qlwice2 : c c C The net heat flux into the ocean is given by: c c Qnet = rlh+sh+qlw+qsw+qios c C The net fresh water flux into the ocean is given by: c c Fnet = fac*rlh + pp + brne c c whith fac=-1/(Qlat*rho_ocean) c where Qlat: latent heat of fusion [2.5e6 J/kg] c and rho_ocean: reference density for ocean [1000 kg/m^3] c c Ice properties: c The mean ice temperature can be calculated by c ticem=thice/hice*cice+tfreeze c The actual ice thickness can be calculated by c hicea=hice/cice c The ice volume by c icevol=hice*dyd*dxd c c Several parameter are taken from include file !!?? c include 'amlice.h' dimension qisw(mx,my),cldfr(mx,my),wspd(mx,my),u(mx,my),v(mx,my), + ppi(mx,my) dimension lsm(mx,my),dyd(my),dxd(mx,my),q(mx,my),t(mx,my) dimension sst(mx,my),sss(mx,my) dimension rlh(mx,my),sh(mx,my),qlw(mx,my),qsw(mx,my),pp(mx,my) dimension qa(mx,my),th(mx,my),rh(mx,my) dimension hice(mx,my), cice(mx,my), thice(mx,my), tsnw(mx,my) dimension rlhi(mx,my),shi(mx,my),qlwi(mx,my),qswi(mx,my), + brne(mx,my),qios(mx,my) dimension rlc0ice(mx,my),cpc0ice(mx,my), + qlwice1(mx,my),qlwice2(mx,my) c c c Additional arrays internally used only c dimension aiflux(4) c if (mx.lt.nx.or.my.lt.ny) stop 'atmosice: dimens. mx,my to small' c c avoid problems with zero time step c tstep=max(1e-8,tstep) c c First call atmospheric heat fluxes assuming to old ice concentration c call htfluxi(tstep,tmonth,sst,tsnw,cice,u,v,wspd,lsm,q,t,cldfr, + sh,rlh,qlw,qa,th,rh,rlc0ice,cpc0ice,qlwice1,qlwice2, + slat,dxd,dyd,nx,ny,mx,my,pbl_eddy,pbl_eddy_nlat,pbl_eddy_nwidth, + pbl_eddy_slat,pbl_eddy_swidth,pbl_eddy_xbox,pbl_eddy_ybox,c0_fac,iperx) c Loop to call 1D thermodynamic ice model c do 10 ix=1,nx do 10 iy=1,ny c c check for land or no ice possible c if (sst(ix,iy).gt.ssticemax .or. lsm(ix,iy).eq.1) then c c set output for no-ice situation c pp(ix,iy)=ppi(ix,iy) qsw(ix,iy)=qisw(ix,iy)*(1-albedoocean) qios(ix,iy)=0. brne(ix,iy)=0. hice(ix,iy)=0. thice(ix,iy)=0. tsnw(ix,iy)=sst(ix,iy) cice(ix,iy)=0. rlhi(ix,iy)=0. shi(ix,iy)=0. qlwi(ix,iy)=0. qswi(ix,iy)=0. else c c -atmos-ice flux parts (sensible, latent, shortwave, longwave) c full flux given by aiflux(1)+aiflux(2)*tsnow+aiflux(3)*qs(tsnow) c t1=qlwice1(ix,iy)-qlwice2(ix,iy)*th(ix,iy) t2=-rlc0ice(ix,iy)*qa(ix,iy) t3=-cpc0ice(ix,iy)*th(ix,iy) albedo=albedoocean+(albedoice-albedoocean)* + exp((tsnw(ix,iy)-tfreeze)/albedof) aiflux(1)=t1+t2+t3+qisw(ix,iy)*(1-albedo) aiflux(2)=cpc0ice(ix,iy)+qlwice2(ix,iy) aiflux(3)=rlc0ice(ix,iy) c -atmos-ice precip aiflux(4)=ppi(ix,iy) c c call 1D ice-thermodynamic model c call icethermo(tsnw(ix,iy),tstep, 1 sst(ix,iy),sss(ix,iy),aiflux,thice(ix,iy),hice(ix,iy), 2 cice(ix,iy),qios(ix,iy),brne(ix,iy),niter,qsice) c c set atmos-ocean fluxes up c sh(ix,iy)=sh(ix,iy)*(1-cice(ix,iy)) rlh(ix,iy)=rlh(ix,iy)*(1-cice(ix,iy)) qlw(ix,iy)=qlw(ix,iy)*(1-cice(ix,iy)) qsw(ix,iy)=qisw(ix,iy)*(1-cice(ix,iy))*(1-albedoocean) pp(ix,iy)=ppi(ix,iy)*(1-cice(ix,iy)) c c set atmos-ice fluxes up c shi(ix,iy)=cpc0ice(ix,iy)*(tsnw(ix,iy)-th(ix,iy))*cice(ix,iy) rlhi(ix,iy)=rlc0ice(ix,iy)*(qsice-qa(ix,iy))*cice(ix,iy) qlwi(ix,iy)=(qlwice1(ix,iy)+qlwice2(ix,iy)* 1 (tsnw(ix,iy)-th(ix,iy)))*cice(ix,iy) c c compute temperature dependent albedo c albedo decays exponentially from ocean to ice value for snow c temperatures below freezing with a decayscale given by c albedof [1/K] c albedo=albedoocean+(albedoice-albedoocean)* + exp((tsnw(ix,iy)-tfreeze)/albedof) qswi(ix,iy)=qisw(ix,iy)*cice(ix,iy)*(1-albedo) c endif 10 continue return end c c=============================================================================== c c subroutine icethermo(tsnw,tstep, + sst,sss,aiflux,thice,hice,cice,qios,brne,iter,qsice) c c======================================================================= c A simple 1D-thermodynamic ice model. = c Closely follows the physics of Hibler and Oberhuber = c = c Solves for ice temperature, growth and melt. = c A timestep, and atmospheric fluxes are specified, hopefully by an = c atmospheric model. = c = c======================================================================= c = c tair = c = c ---- qas (atm-ice flx) ---------------- = c tsnw |hsnow = c ---- qsi (=qas) ----------------------- = c tice | = c | = c qif (ice) | = c tfreeze | = c -----qio (ice-ocean)-----------------------------------------------= c SST = c======================================================================= c final ice-ocean heat flux : qios [W/m**2)] = c final freshwater flux : brne [m/s] = c grid mean ice thickness : hice [m] = c fraction of grid ice covered : cice fraction = c mean 'heat' content if ice : thice [K*m] = c = c note that hsnow is fixed and has no heat content .... = c = c====================================================================== c --- Set some constants for the run from include file !!?? include 'amlice.h' dimension aiflux(4) c aiflux(1-3) is used ti evaluate the ice-atmos heat flux c aiflux(4) is precip c dqmax is the maximum missfit in W/m^2 between ice-snow and atmos-snow flux c====================================================================== c --- ice-ocean heat flux c modeled as a conductive heat flux c tkocean :ocean flux coeff [W/m^2K] could be rho(0)*Cp*hmix/dt c but is assumed to be constant (typical value ~ 5000 w/m^2K) c tfreeze : freezing temeprature assumed to be fixed (-1.8 C) qio=tkocean*(sst-tfreeze) c --- inititalize some variables dqf=100. iter=0 c== check if ice exist already if (hice.gt.hicemin) then c --- get old ice temperature c hice: represents the 'heat' content in units [Cm] c tice: temperature at top of ice sheet. c Assuming the ice bottom is at tfreeze, and linear temp. profile tice=2*thice/hice*cice+tfreeze c --- qsi is the heat flux at the snow ice interface qsi=(tksnow/hsnow)*(tice-tsnw) c qif is the heat flux through the ice qif=(tkice/hice*cice)*(tfreeze-tice) c tkice is the ice flux coeff in W/m^2K c --- start iteration loop tconv=tconvin do while ((abs(dqf).gt.dqmax).and.(iter.lt.itermax)) iter=iter+1 c c tconv is convergence factor for the higly nonlinear flux coupling c tconv=tconv*tconvgr tconv=min(tconv,tconvmax) c --- atmos-snow heat flux c here we need to know what is done in the atmospheric model c and do the same thing each iteration to get the propoer fluxes. c First solve for the saturation humidity over ice qsice=0.622*6.11/1000*exp(17.67*(tsnw-273.15)/ + (tsnw-273.15+243.5)) qsice=qsice*10.**(0.00422*(tsnw-273.15)) c --- get current atmos-snow flux qas=aiflux(1)+aiflux(2)*tsnw+aiflux(3)*qsice c --- if divergent try to reduce convergence factor if (abs(dqf).lt.abs(qas-qsi)) then tconv=tconv*0.2 endif c --- get missmatch between atmos-snow and snow-ice flux dqf=qas-qsi c --- relax snow-icetop heat flux toward atmos-snow heat flux qsi=qsi+tconv*dqf c --- solve for appropriate snow temperature tsnw=tice-qsi*hsnow/tksnow c --- if snow temperature is at freezing ignore snow layer if (tsnw-tfreeze.gt.0.5) then dqf=0. tsnw=tfreeze qsice=0.622*6.11/1000*exp(17.67*(tsnw-273.15)/ 1 (tsnw-273.15+243.5)) qsice=qsice*10.**(0.00422*(tsnw-273.15)) qas=aiflux(1)+aiflux(2)*tsnw+aiflux(3)*qsice qsi=qas endif enddo c end of iteration loop c c change heat content of ice due to conductive fluxes c thice=thice-(tstep/(cpice*rhoice))*(qsi-qif)*cice c c check for melting due to conductive fluxes c if (thice.gt.0) then qmelt=thice * cpice * rhoice / tstep thice=0. else qmelt = 0. endif else c c --- when no ice existed do the following c tsnw=tfreeze if(qio.lt.0) then cice=0.3 else cice=0.0 endif qsi=0.0 qif=0.0 tice=0. endif c ===== find out how much ice will be grown or melted at this timestep === c --- get ice growth/melt at top due to precip and surface melt ppi=cice*aiflux(4)*rhowater/rhoice dhpdt=ppi-qmelt/(rhoice*hfusionice) c c get change in ice thickness at bottom of ice dhdt=((qif-qio)*cice)/(rhoice*hfusionice) +dhpdt c dont melt more ice than existed dhdt=max(dhdt,-hice/tstep) dhpdt=max(dhpdt,-hice/tstep) c save new ice thickness hice=hice+dhdt*tstep c --- get (relatively) fresh water flux, brne = -(dhdt-ppi)*rhoice/rhowater*(1-sice/sss) c --- save ice-ocean heat flux qio*cice backed from ice growth qios = (dhpdt-dhdt)*rhoice*hfusionice + qif*cice c c --- no heatcontent for ice thinner than hicemin c if (hice.lt.hicemin) then thice=0.0 qif=0.0 qsi=0.0 qio=0.0 tsnw=tfreeze endif c ======== check if ice has grown or melted and change concentration c --- forget about the last millimeter of ice if (hice.lt.1.e-3) then hice=0.0 cice=0.0 else c --- change concentration; dq is the concentration change. if (dhdt.gt.0) then c --- freezing dq=dhdt*tstep/hq*(1-cice) else c --- melting dq=dhdt*tstep/(hf*hice)*cice endif c --- change ice concentration and limit to bounds cice=min(cice+dq,cicemax) cice=max(cice,0.1) endif return end c ====================================================================== c LIST OF VARIABLES: c ----------------- c thice :ice temp[C m] c hice :ice thickness [m] c cice :ice concentration [0-1] c rhoice :ice density 910 [kg/m^3] c tksnow :snow conductivity 0.33 [W/mK] c tkice :ice conductivity 2.0 [W/mK] c hsnow :snow thickness [m] c hicemin :minimum ice thickness [m] c hfusionic :latent heat of fusion for ice 3.34e5 [J/kg] c tfreeze :freezing temperature -1.8 [C] c hq :convert dh->dq (thickness to extent) freezing 0.25 [m] c hf :convert dh->dq (thickness to extent) melting 2.0 c cpice :specific heat 2090 [J/kgK] c tkocean :ocean flux coeff [W/M^2K] could be rho0*cp*hmix/dt c tsnw :snow temperature c sice :ice salinity c itermax :maximum nuber of ice iterations c qio :ice-ocean heat flux. c qif :icetop-icebottom heat flux. c qsi :snow-icetop heat flux. c qas :atmos-snow heat flux. c qios :net ice-ocean flux (no lead contribution) c brne :net ice-ocean freshwater flux c============================================================================ c c subroutine htfluxi (tstep,tmonth,sst,tice,fice,u,v,wspd,lsm,q,t,cldfr, $ sh,rlh,qlw,qa1,th,rh,rlc0ice,cpc0ice,qlwice1,qlwice2, $ slat,dxd,dyd,nx,ny,mx,my,coef_eddy,offn,stdn,offs,stds, $ xbox,ybox,c0_fac,iperx) 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 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 Also, it should be noted that the code is set up for the c case of a basin bounded at the east and west. It hence c cannot deal with the part of the Southern Ocean that goes c through the Drake passage, or the part of the Arctic north c of Greenland. In both case the matrices would no longer be c tridiagonal. We will change this but we're talking months c and months (years?) here. c c The inputs are: c c All temperatures in K c c tstep = time step (sec) c tmonth = month of the year (sort of) c sst = array containing the model or observed SST c tice = array containing the model or observed sea ice temperature c fice = array containing the model or observed sea ice fraction 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 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 <= mx c ny = number of y grid points <= my c mx = x grid dimension c my = y grid dimension c c The outputs are: c c sh = array containing the sensible heat flux for water(W/m^2) c rlh = array containing the latent heat flux for water(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 for water c rh = relative humidity as a fraction c rlc0ice = rl*c0ice(i,j)*wspd(i,j)*rhoa c cpc0ice = cp*c0ice(i,j)*wspd(i,j)*rhoa c qlwice1 = factors*th(i,j)**4 c qlwice2 = factors*th(i,j)**3 c c All fluxes are for the ocean fraction only. An ice model will c calculate the fluxes over sea ice using the air temperature and c air humidity derived here. c c The longwave flux over ice can be got from: c c qlwice=qlwice1+qlwice2*(tice-th(i,j)) 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 c constants: rl=latent heat of water. rlice . . . ice c cp = specific heat of water at constant pressure c r=univ. gas constant, stef=stefan bolz.'s const. implicit real*4(a-h,o-z),integer(i-n) c make sure parameter nxx, nyx are bigger or equal to mx,my c parameter (nxx=??,nyx=??) parameter (nxx=300,nyx=300) c!!! don't change nxx,nyx without changing maxdim in ADVDIFQ1DX and ADVDIFQ1D dimension sst(mx,my),u(mx,my),v(mx,my),wspd(mx,my),q(mx,my), $ t(mx,my),rlh(mx,my),sh(mx,my),lsm(mx,my),qa1(mx,my), $ th(mx,my),qlw(mx,my),cldfr(mx,my),dyd(my), $ dxd(mx,my),rh(mx,my),tice(mx,my),fice(mx,my), $ rlc0ice(mx,my),cpc0ice(mx,my),qlwice1(mx,my),qlwice2(mx,my) dimension up(nxx,nyx),vp(nxx,nyx),thv(nxx,nyx), $ thve(nxx,nyx),qs(nxx,nyx),dy(nyx), $ thvs(nxx,nyx),pnuxp(nxx,nyx), $ pnuyp(nxx,nyx),c0(nxx,nyx),qe(nxx,nyx), $ dx(nxx,nyx),qa(nxx,nyx), $ c0thv(nxx,nyx),c0q(nxx,nyx),qsice(nxx,nyx), $ c00(nxx,nyx),c0ice(nxx,nyx),thvst(nxx,nyx) $ ,dtdx(nxx,nyx),dtdy(nxx,nyx) c# $ ,dtdx_ave(nxx,nyx),dtdy_ave(nxx,nyx) integer idim(2) logical advec if(nxx.lt.nx .or. nyx.lt.ny) then write(*,*) 'arrays in subroutine htfluxi are dimensioned less $ than ny and nx set in calling routine' write(*,*) 'nxx,nx=',nxx,nx write(*,*) 'nyx,ny=',nyx,ny stop endif c advec=.true. implements advection advec=.true. jstart=2 jend=ny-1 c set model parameters pnu=0.4e+7 c# pnu=0. delta=.25 c#try delta=.18 pml=6000. depth=600. betav=0.17 qrad=-2./86400. c set constants r=287.04 psfc=100000. rl=2.5e+6 rlice=2.834e+6 cp=1004. rhoa=1.225 stef=5.6696e-8 eps=0.97 idim(1)=nx idim(2)=ny c determine grid spacing in m conv=2.*3.14/360. radius=6.37e+6 dy(1)=radius*dyd(1)*conv rlat=slat*conv do 1 i=1,nx dx(i,1)=conv*radius*cos(rlat)*dxd(i,1) 1 continue do 2 j=2,ny dy(j)=conv*radius*dyd(j) rlat=rlat+dyd(j)*conv do 2 i=1,nx dx(i,j)=conv*radius*cos(rlat)*dxd(i,j) 2 continue 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. fac=.622*6.11/1000. do 24 j=1,ny do 24 i=1,nx c00(i,j)=0.0014*c0_fac c0ice(i,j)=0.0028*c0_fac if(lsm(i,j) .eq. 0) then qs(i,j)=fac*exp(17.67*(sst(i,j)- $ 273.15)/(sst(i,j)-273.15+243.5)) if(fice(i,j).gt.0.) then qsice1=fac*exp(17.67*(tice(i,j)- $ 273.15)/(tice(i,j)-273.15+243.5)) qsice(i,j)=qsice1*10.**(0.00422*(tice(i,j)-273.15)) endif endif 24 continue do j = 2, ny-1 do i = 2, nx-1 dtdx(i,j) = (th(i,j) - th(i-1,j))/dx(i,j) enddo enddo do j = 2, ny-1 do i = 2, nx-1 dtdy(i,j) = (th(i,j) - th(i,j-1))/dy(j-1) enddo enddo c Richard says the coef_eddy should be 0.18e5 [m^2/s/K] eddy_coefn = coef_eddy * (1. + cos((tmonth-1.)*3.14159/6.))/2. eddy_coefs = coef_eddy * (1. - cos((tmonth-1.)*3.14159/6.))/2. rlat=slat do j = 1, ny c# eddy_coef = eddy_coefn c# if (rlat.lt.0) eddy_coef = eddy_coefs ygn=(rlat-offn)/stdn ygs=(rlat-offs)/stds eddy_coef = eddy_coefn* exp(- (ygn**2) / 2.) + * eddy_coefs* exp(- (ygs**2) / 2.) c# eddy_coef2 = eddy_coef * (1. - cos (4.*abs(rlat*conv)) ) c# write(98,*) rlat,eddy_coef rlat = rlat + dyd(j) do i = 1, nx accumx = 0. accumy = 0. xsum = 0. im = min(nx-1,i) num = 0 do while (xsum .lt. xbox) jp = j jm = j ysum = 0. do while (ysum .lt. ybox) ysum = ysum + dyd(jp) + dyd(jm) accumx = accumx + dtdx(im,jp) accumy = accumy + dtdy(im,jp) accumx = accumx + dtdx(im,jm) accumy = accumy + dtdy(im,jm) jm = max(2,jm-1) jp = min(ny-1,jp+1) num = num + 2 enddo im = max(2,im-1) xsum = xsum + dxd(im,jm) enddo c# dtdx_ave(i,j) = eddy_coef* accumx/ num c# dtdy_ave(i,j) = eddy_coef* accumy/ num u(i,j) = u(i,j) + eddy_coef* accumx/ num v(i,j) = v(i,j) + eddy_coef* accumy/ num enddo enddo pfac=psfc/(psfc-.5*pml) c#nhn iter=1 c#nhn itermx=3 c#nhn 99 do 25 j=1,ny do 25 j=1,ny do 25 i=1,nx c#nhn if(iter.gt.1 .and. (thv(i,j).gt.thvst(i,j))) then c#nhn c00(i,j)=.00075*c0_fac c#nhn endif if(lsm(i,j).eq.1) then thve(i,j)=t(i,j)*(1.+.61*q(i,j)) qe(i,j)=q(i,j) th(i,j)=t(i,j) qa1(i,j)=q(i,j) else w0=wspd(i,j)*pml/depth thvs(i,j)=sst(i,j)*(1.+.61*qs(i,j)) if(fice(i,j).eq.0.) then thve(i,j)=thvs(i,j)+pml*qrad/((1.+betav)*c00(i,j)*w0) qe(i,j)=qs(i,j)/(1.+delta) thvst(i,j)=thvs(i,j) else thvice=tice(i,j)*(1.+.61*qsice(i,j)) c0thv(i,j)=fice(i,j)*c0ice(i,j)+ $ (1.-fice(i,j))*c00(i,j) c0q(i,j)=fice(i,j)*c0ice(i,j)+ $ (1.-fice(i,j))*c00(i,j) thvst(i,j)=(fice(i,j)*c0ice(i,j)*thvice+ $ (1.-fice(i,j))*c00(i,j)*thvs(i,j))/c0thv(i,j) qst=(fice(i,j)*c0ice(i,j)*qsice(i,j)+ $ (1.-fice(i,j))*c00(i,j)*qs(i,j))/c0q(i,j) thve(i,j)=thvst(i,j)+pml*qrad/((1.+betav)*c0thv(i,j)*w0) qe(i,j)=qst/(1.+delta) endif endif 25 continue c Set equilibrium values to observed at northernmost and southernmost c points if they are open ocean. This is required because c 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. c ttoth is a conversion from observed 1000mb temperature to surface c temperature using a dry adiabat. This is here 'cos my input data c was for 1000mb temperature not surface temperature and 'cos the c observed SLP beyond the extremes of the grid was greater than c 1000mb. You can do what you want but be careful! ttoth=(1017./(1000.))**(r/cp) do 26 i=1,nx do 27 j=1,jstart-1 if(lsm(i,j).eq.0) then 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)*ttoth*(1.+.61*q(i,j)) qa(i,j)=q(i,j) thv(i,j)=thve(i,j) th(i,j)=t(i,j)*ttoth endif 27 continue do 26 j=jend+1,ny if(lsm(i,j).eq.0) then 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)*ttoth*(1.+.61*q(i,j)) qa(i,j)=q(i,j) thv(i,j)=thve(i,j) th(i,j)=t(i,j)*ttoth endif 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 w0=wspd(i,j)*pml/depth if(lsm(i,j).eq.1) then up(i,j)=0. vp(i,j)=0. pnuxp(i,j)=0. pnuyp(i,j)=0. else if(fice(i,j).eq.0.) then c0(i,j)=c00(i,j) else c0(i,j)=c0thv(i,j) endif 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).eq.1 .or. lsm(im1,j).eq.1 $ .or. lsm(i,jp1).eq.1 .or. lsm(i,jm1).eq.1 $ .or. lsm(ip2,j).eq.1 .or. lsm(im2,j).eq.1 $ .or. lsm(i,jp2).eq.1 .or. lsm(i,jm2).eq.1 ) 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(i,j).gt.0.) then i1=i else i1=i+1 if(i.eq.nx) i1=nx endif up(i,j)=u(i,j)*pml/((1.+betav)*c0(i,j)*w0*dx(i1,j)) if(v(i,j).gt.0.) then vp(i,j)=v(i,j)*pml/((1.+betav)*c0(i,j)*w0*dy(j-1)) else vp(i,j)=v(i,j)*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 adv2Deq1mp(iperx,idim,nxx,nyx,up,vp,pnuxp,pnuyp,thve,thv) c#new call advdiff(iperx,idim,nxx,nyx,up,vp,pnuxp,pnuyp,thve,thv) c repeat one time c#nhn iter=iter+1 c#nhn 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 w0=wspd(i,j)*pml/depth if(lsm(i,j).eq.0) then if(fice(i,j).eq.0.) then c0(i,j)=c00(i,j) else c0(i,j)=c0q(i,j) endif 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).eq.1 .or. lsm(im1,j).eq.1 $ .or. lsm(i,jp1).eq.1 .or. lsm(i,jm1).eq.1 $ .or. lsm(ip2,j).eq.1 .or. lsm(im2,j).eq.1 $ .or. lsm(i,jp2).eq.1 .or. lsm(i,jm2).eq.1 ) 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(i,j).gt.0.) then i1=i else i1=i+1 if(i.eq.nx) i1=nx endif up(i,j)=u(i,j)*pml/((1.+delta)*c0(i,j)*w0*dx(i1,j)) if(v(i,j).gt.0.) then vp(i,j)=v(i,j)*pml/((1.+delta)*c0(i,j)*w0*dy(j-1)) else vp(i,j)=v(i,j)*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 adv2Deq1mp(iperx,idim,nxx,nyx,up,vp,pnuxp,pnuyp,qe,qa) c calculate theta from theta_V and q c calculate fluxes of sensible and latent heat do 30 j=1,ny do 30 i=1,nx if(lsm(i,j).eq. 0) then rlh(i,j)=rhoa*rl*c00(i,j)*wspd(i,j)*(qs(i,j)-qa(i,j)) th(i,j)=thv(i,j)/(1.+.61*qa(i,j)) sh(i,j)=rhoa*cp*c00(i,j)*wspd(i,j)*(sst(i,j)-th(i,j)) c_new stuff follows sstk = sst(i,j) if ( sstk .gt. 301.15 ) then aout =.4 else aout =.8 endif qlw(i,j)=eps*stef*sstk**3*(sstk*(.417-.0486* $ sqrt(abs(qa(i,j))*1000./.622))* $ (1.-aout*cldfr(i,j)**2) + 4.*(sstk-th(i,j))) c_new stuff ends c_orig qlw(i,j)=eps*stef*(th(i,j)**4.)*(.39-.05* c_orig $ sqrt(abs(qa(i,j))*1000./.622)) c_orig $ *(1.-.55*cldfr(i,j)) + 4.*eps*stef* c_orig $ (th(i,j)**3.)*(sst(i,j)-th(i,j)) rlc0ice(i,j)=rhoa*rlice*c0ice(i,j)*wspd(i,j) cpc0ice(i,j)=rhoa*cp*c0ice(i,j)*wspd(i,j) qlwice1(i,j)=eps*stef*(th(i,j)**4.)*(.39-.05* $ sqrt(abs(qa(i,j))*1000./.622)) $ *(1.-.55*cldfr(i,j)) qlwice2(i,j)=4.*eps*stef*(th(i,j)**3.) qa1(i,j)=qa(i,j) endif qsatair=fac*exp(17.67*(th(i,j)-273.15)/(th(i,j)-273.15+243.5)) rh(i,j)=qa(i,j)/qsatair !relative humidity 30 continue return end c============================================================================ c SUBROUTINE adv2Deq1m(IDIM,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) INTEGER IDIM(2) C variables are dimensioned with X first MX = IDIM(1) MY = IDIM(2) NXSKP = 1 NYSKP = NX C does X advection C loops over all latitudes IX = 1 DO IY = 1 , MY CALL ADVDIFQ1DX(UP(1,IY),NUXP(1,IY),MX,QE(1,IY),QE(1,IY), * QE(MX,IY),QA(1,IY)) END DO C does Y advection C loops over all longitudes IY = 1 DO IX = 1 , MX C boundary conditions QLEFT = QE(IX,1) QRIGHT = QE(IX,MY) CALL ADVDIFQ1D(VP(IX,1),NUYP(IX,1),MY,QA(IX,1), * QLEFT,QRIGHT,NYSKP,QA(IX,1)) END DO RETURN END c============================================================================ subroutine adv2Deq1mp (iper,IDIM,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) INTEGER IDIM(2) C variables are dimensioned with X first MX = IDIM(1) MY = IDIM(2) NXSKP = 1 NYSKP = NX C does Y advection C loops over all longitudes IY = 1 DO IX = 1 , MX C boundary conditions QLEFT = QE(IX,1) QRIGHT = QE(IX,MY) CALL ADVDIFQ1D(VP(IX,1),NUYP(IX,1),MY,QE(IX,1), * QLEFT,QRIGHT,NYSKP,QA(IX,1)) END DO C does X advection C loops over all latitudes IX = 1 if (iper.eq.0) then DO IY = 1 , MY CALL ADVDIFQ1DX(UP(1,IY),NUXP(1,IY),MX,QA(1,IY),QA(1,IY), 1 QA(MX,IY),QA(1,IY)) END DO else DO IY = 1 , MY CALL ADVDIFQ1DXP(UP(1,IY),NUXP(1,IY),MX,QA(1,IY),QA(1,IY), 1 QA(MX,IY),QA(1,IY)) END DO endif RETURN END c============================================================================ SUBROUTINE ADVDIFQ1DXP(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)+NU2(1))*QRIGHT BC(1) = 1. + U2(1) + 2*NU2(1) CC(1) = - NU2(1) ELSE C boundary condition doesn't matter QE(1) = QE0(1)+ NU2(1)*QRIGHT BC(1) = 1. - U2(1) + 2*NU2(1) CC(1) = U2(1) - NU2(1) ENDIF C does right point IF(U2(NX).GT.0.0)THEN C boundary condition doesn't matter QE(NX) = QE0(NX) + NU2(NX)*QLEFT BC(NX) = 1. + U2(NX) + 2*NU2(NX) AC(NX) = -U2(NX) - NU2(NX) ELSE C boundary condition matters QE(NX) = QE0(NX) - (U2(NX)-NU2(NX))*QLEFT BC(NX) = 1.0 - U2(NX) + 2*NU2(NX) AC(NX) = - NU2(NX) ENDIF CALL TRIDAGQA(AC,BC,CC,QE,CC,QE,QA,1,NX) RETURN END