subroutine MFLUX2 (fxh, fxe, fxmx, fxmy, cdm, cdm2,xxfh) !------------------------------------------------------------------------ ! ! MFLUX2 computes surface fluxes of momentum, heat,and moisture ! using monin-obukhov. the roughness length "z0" is prescribed ! over land and over ocean "z0" is computed using charnocks formula. ! the universal functions (from similarity theory approach) are ! those of hicks ! !------------------------------------------------------------------------ use allocate_mod include 'RESOLUTION.h' include 'PARAMETERS.h' include 'STDUNITS.h' include 'FLAGS.h' include 'BKINFO.h' include 'CONSLEV.h' include 'CONMLEV.h' include 'ESTAB.h' include 'FILEC.h' include 'FILEPC.h' include 'GDINFO.h' include 'IIWD.h' include 'LIMIT.h' include 'NUM.h' include 'QLOGS.h' include 'STORE.h' include 'TIME.h' include 'WINDD.h' include 'ZLDATA.h' common /dif/ ihsflp common /idiag/ idiagn !----------------------------------------------------------------------- ! user interface variables !----------------------------------------------------------------------- real, intent (out), dimension (ibeg:iend) :: fxh real, intent (out), dimension (ibeg:iend) :: fxe real, intent (out), dimension (ibeg:iend) :: fxmx real, intent (out), dimension (ibeg:iend) :: fxmy real, intent (out), dimension (ibeg:iend) :: cdm real, intent (out), dimension (ibeg:iend) :: cdm2 !----------------------------------------------------------------------- ! internal variables !----------------------------------------------------------------------- parameter (icntx = 30) data amask/ -98.0/ integer, dimension(imx) :: ifz integer, dimension(imx) :: indx integer, dimension(imx) :: istb integer, dimension(imx) :: it integer, dimension(imx) :: iutb real, dimension(imx) :: aap real, dimension(imx) :: bq1 real, dimension(imx) :: bq1p real, dimension(imx) :: delsrad real, dimension(imx) :: ecof real, dimension(imx) :: ecofp real, dimension(imx) :: estso real, dimension(imx) :: estsop real, dimension(imx) :: fmz1 real, dimension(imx) :: fmz10 real, dimension(imx) :: fmzo1 real, dimension(imx) :: fhzo1 real, dimension(imx) :: foft real, dimension(imx) :: foftm real, dimension(imx) :: frac real, dimension(imx) :: land real, dimension(imx) :: pssp real, dimension(imx) :: qf real, dimension(imx) :: rdiff real, dimension(imx) :: rho real, dimension(imx) :: rib real, dimension(imx) :: rkmaxp real, dimension(imx) :: rstso real, dimension(imx) :: rstsop real, dimension(imx) :: sf10 real, dimension(imx) :: sfm real, dimension(imx) :: sfh real, dimension(imx) :: sfzo real, dimension(imx) :: sgzm real, dimension(imx) :: slwa real, dimension(imx) :: szeta real, dimension(imx) :: szetam real, dimension(imx) :: t1 real, dimension(imx) :: t2 real, dimension(imx) :: tab1 real, dimension(imx) :: tab2 real, dimension(imx) :: tempa1 real, dimension(imx) :: tempa2 real, dimension(imx) :: theta real, dimension(imx) :: thetap real, dimension(imx) :: tsg real, dimension(imx) :: tsm real, dimension(imx) :: tsp real, dimension(imx) :: tss real, dimension(imx) :: ucom real, dimension(imx) :: uf10 real, dimension(imx) :: ufh real, dimension(imx) :: ufm real, dimension(imx) :: uft real, dimension(imx) :: ufzo real, dimension(imx) :: ugzm real, dimension(imx) :: uzeta real, dimension(imx) :: uzetam real, dimension(imx) :: vcom real, dimension(imx) :: vrtkx real, dimension(imx) :: vrts real, dimension(imx) :: wind real, dimension(imx) :: windp real, dimension(imx) :: xxfh real, dimension(imx) :: xxfm real, dimension(imx) :: xxfm2 real, dimension(imx) :: xxsh real, dimension(imx) :: z10 real, dimension(imx) :: zeta real, dimension(imx) :: zkmax real, dimension(imx) :: pss real, dimension(imx) :: tstar real, dimension(imx) :: ukmax real, dimension(imx) :: vkmax real, dimension(imx) :: tkmax real, dimension(imx) :: rkmax real :: yz,y1,y2,y3,y4,windmks,znotm character*10 routine routine = 'mflux2' yz= 0.0001344 y1= 3.015e-05 y2= 1.517e-06 y3= -3.567e-08 y4= 2.046e-10 !------------------------------------------------------------------------ ! set water availability constant "ecof" and land mask "land". ! limit minimum wind speed to 100 cm/s !------------------------------------------------------------------------ j = IFIX(tjloc(i1)) do i = 1,imx z10(i) = 1000. pss(i) = pspc(i) tstar(i) = tstrc(i) ukmax(i) = upc(kmax,i) vkmax(i) = vpc(kmax,i) tkmax(i) = tpc(kmax,i) rkmax(i) = rpc(kmax,i,1) enddo do i = i1,i2 windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i)) wind (i) = amax1(windp(i),100.) if (zoc(i) .GT. 0.0) then ecof(i) = wetc(i) land(i) = 1.0 zot (i) = zoc(i) else windmks=wind(i)*.01 ecof(i) = wetc(i) land(i) = 0.0 if(wind(i) .lt. 1250.0 ) then zoc(i)=-(0.0185/9.8*(7.59e-8*wind(i)**2+ & 2.46e-4*wind(i))**2)*100. endif if(wind(i) .ge. 1250.0 .and. wind(i) .lt. 3000.0) & zoc(i)=-(.000739793 * wind(i) -0.58)/10 if(wind(i) .ge. 3000.0) then znotm=yz+windmks*y1+windmks**2*y2+windmks**3*y3+ & windmks**4*y4 !powell 2003 zoc(i) = -100.*znotm endif if(windmks .le. 7.0) then znott = (0.0185/9.8*(7.59e-8*wind(i)**2+ & 2.46e-4*wind(i))**2) zot(i) = -100* znott endif if(windmks .gt. 7.0 .and. windmks .lt. 20.0 ) then cccc hwrf if(windmks .gt. 7.0) then znott=0.2375*exp(-0.5250*windmks)+ & 0.0025*exp(-0.0211*windmks) zot(i) = -znott endif if(windmks .ge. 20.0 .and. windmks .lt. 40.0 ) then znott=1.634 * (10.0**(-14)) * (windmks ** 6.8905) zot(i) = -100.0 * znott endif if(windmks.ge.40.) zot(i) = zoc(i) endif c if(nstep.eq.0) zot(i) = zoc(i) if (zoc(i) .LT. amask) & zoc(i) = -0.0185*0.001*wind(i)*wind(i)*og if (zot(i) .LT. amask) & zot(i) = -0.0185*0.001*wind(i)*wind(i)*og !------------------------------------------------------------------------ ! where necessary modify zo values over ocean. !------------------------------------------------------------------------ enddo !------------------------------------------------------------------------ ! define constants: ! a and c = constants used in evaluating universal function for ! stable case ! ca = karmen constant ! cm01 = constant part of vertical integral of universal ! function; stable case ( 0.5 < zeta < or = 10.0) ! cm02 = constant part of vertical integral of universal ! function; stable case ( zeta > 10.0) !------------------------------------------------------------------------ en = 2. c = .76 a = 5. ca = .4 cmo1 = .5*a - 1.648 cmo2 = 17.193 + .5*a - 10.*c boycon = .61 do i = i1,i2 theta(i) = tkmax(i)*rqc9 vrtkx(i) = 1.0 + boycon*rkmax(i) zkmax(i) = rgas*tkmax(i)*qqlog(kmax)*og enddo !------------------------------------------------------------------------ ! get saturation mixing ratios at surface !------------------------------------------------------------------------ do i = i1,i2 tsg (i) = tstar(i) tab1 (i) = tstar(i) - 153.16 it (i) = IFIX(tab1(i)) tab2 (i) = tab1(i) - FLOAT(it(i)) t1 (i) = tab(it(i) + 1) t2 (i) = table(it(i) + 1) estso(i) = t1(i) + tab2(i)*t2(i) psps1 = (pss(i) - estso(i)) if(psps1 .EQ. 0.0)then psps1 = .1 endif rstso(i) = 0.622*estso(i)/psps1 vrts (i) = 1. + boycon*ecof(i)*rstso(i) enddo !------------------------------------------------------------------------ ! check if consideration of virtual temperature changes stability. ! if so, set "dthetav" to near neutral value (1.0e-4). also check ! for very small lapse rates; if ABS(tempa1) <1.0e-4 then ! tempa1=1.0e-4 !------------------------------------------------------------------------ do i = i1,i2 tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i) tempa2(i) = tempa1(i)*(theta(i) - tstar(i)) if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4 tab1(i) = ABS(tempa1(i)) if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4 !------------------------------------------------------------------------ ! compute bulk richardson number "rib" at each point. if "rib" ! exceeds 95% of critical richardson number "tab1" then "rib = tab1" !------------------------------------------------------------------------ rib (i) = ggg*zkmax(i)*tempa1(i)/ & (tkmax(i)*vrtkx(i)*wind(i)*wind(i)) tab2(i) = ABS(zoc(i)) tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i))) if (rib(i) .GT. tab1(i)) rib(i) = tab1(i) enddo do i = i1,i2 zeta(i) = ca*rib(i)/0.03 enddo !------------------------------------------------------------------------ ! begin looping through points on line, solving wegsteins iteration ! for zeta at each point, and using hicks functions !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! set initial guess of zeta=non - dimensional height "szeta" for ! stable points !------------------------------------------------------------------------ rca = 1./ca enrca = en*rca zog = .0185*og !------------------------------------------------------------------------ ! stable points !------------------------------------------------------------------------ ip = 0 do i = i1,i2 if (zeta(i) .GE. 0.0) then ip = ip + 1 istb(ip) = i endif enddo if (ip .EQ. 0) go to 170 do i = 1,ip szetam(i) = 1.0e+50 sgzm(i) = 0.0e+00 szeta(i) = zeta(istb(i)) ifz(i) = 1 enddo !------------------------------------------------------------------------ ! begin wegstein iteration for "zeta" at stable points using ! hicks(1976) !------------------------------------------------------------------------ do icnt = 1,icntx do i = 1,ip if (ifz(i) .EQ. 0) go to 80 zal1g = ALOG(szeta(i)) if (szeta(i) .LE. 0.5) then fmz1(i) = (zal1g + a*szeta(i))*rca else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then rzeta1 = 1./szeta(i) fmz1(i) = (8.*zal1g + 4.25*rzeta1 - & 0.5*rzeta1*rzeta1 + cmo1)*rca else if (szeta(i) .GT. 10.) then fmz1(i) = (c*szeta(i) + cmo2)*rca endif szetao = ABS(zoc(istb(i)))/zkmax(istb(i))*szeta(i) zal2g = ALOG(szetao) if (szetao .LE. 0.5) then fmzo1(i) = (zal2g + a*szetao)*rca sfzo (i) = 1. + a*szetao else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then rzeta2 = 1./szetao fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - & 0.5*rzeta2*rzeta2 + cmo1)*rca sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2 else if (szetao .GT. 10.) then fmzo1(i) = (c*szetao + cmo2)*rca sfzo (i) = c*szetao endif ! compute heat & moisture parts of zot.. for calculation of sfh szetho = ABS(zot(istb(i)))/zkmax(istb(i))*szeta(i) zal2gh = ALOG(szetho) if (szetho .LE. 0.5) then fhzo1(i) = (zal2gh + a*szetho)*rca sfzo (i) = 1. + a*szetho else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then rzeta2 = 1./szetho fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - & 0.5*rzeta2*rzeta2 + cmo1)*rca sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2 else if (szetho .GT. 10.) then fhzo1(i) = (c*szetho + cmo2)*rca sfzo (i) = c*szetho endif !------------------------------------------------------------------------ ! compute universal function at 10 meters for diagnostic purposes !------------------------------------------------------------------------ if (ngd .EQ. nNEST) then szet10 = ABS(z10(istb(i)))/zkmax(istb(i))*szeta(i) zal10 = ALOG(szet10) if (szet10 .LE. 0.5) then fmz10(i) = (zal10 + a*szet10)*rca else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then rzeta2 = 1./szet10 fmz10(i) = (8.*zal10 + 4.25*rzeta2 - & 0.5*rzeta2*rzeta2 + cmo1)*rca else if (szet10 .GT. 10.) then fmz10(i) = (c*szet10 + cmo2)*rca endif sf10(i) = fmz10(i) - fmzo1(i) endif sfm(i) = fmz1(i) - fmzo1(i) sfh(i) = fmz1(i) - fhzo1(i) sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ & (sfh(i) + enrca*sfzo(i)) fmz = (sgz - szeta(i))/szeta(i) fmzo = ABS(fmz) if (fmzo .GE. 5.0e-5) then sq = (sgz - sgzm(i))/(szeta(i) - szetam(i)) if(sq .EQ. 1) then print*,'NCO: !!! ERROR !!! DIVIDE BY ZERO IN MFLUX2' & ,' (STABLE CASE)' print*,'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i) endif szetam(i) = szeta(i) szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq) sgzm (i) = sgz else ifz(i) = 0 endif 80 continue enddo enddo do i = 1,ip if (ifz(i) .GE. 1) go to 110 enddo ! go to 130 110 continue write(stdout,120) 120 format(2X, ' NON-CONVERGENCE FOR STABLE ZETA IN ROW ') call MPI_CLOSE(1,routine) !------------------------------------------------------------------------ ! update "zo" for ocean points. "zo"cannot be updated within the ! wegsteins iteration as the scheme (for the near neutral case) ! can become unstable !------------------------------------------------------------------------ 130 continue do i = 1,ip szo = zoc(istb(i)) if (szo .LT. 0.0) then wndm=wind(istb(i))*0.01 ckg=0.0185*og szo = - ckg*wind(istb(i))*wind(istb(i))/ & (sfm(i)*sfm(i)) cons_p000001 = .000001 cons_7 = 7. vis = 1.4E-1 ustar = sqrt( -szo / zog) restar = -ustar * szo / vis restar = max(restar,cons_p000001) ! Rat taken from Zeng, Zhao and Dickinson 1997 rat = 2.67 * restar ** .25 - 2.57 rat = min(rat ,cons_7) !constant rat=0. ! zot(istb(i)) = szo * exp(-rat) else zot(istb(i)) = zoc(istb(i)) endif !!cbt zoc(istb(i)) = szo enddo do i = 1,ip xxfm(istb(i)) = sfm(i) c xxfm2(istb(i)) = sfm(i) xxfm2(istb(i)) = sfh(i) xxfh(istb(i)) = sfh(i) xxsh(istb(i)) = sfzo(i) enddo !------------------------------------------------------------------------ ! obtain wind at 10 meters for diagnostic purposes !------------------------------------------------------------------------ if (ngd .EQ. nNEST) then do i = 1,ip wind10(istb(i),j) = sf10(i)*wind(istb(i))/sfm(i) enddo endif !------------------------------------------------------------------------ ! unstable points !------------------------------------------------------------------------ 170 continue iq = 0 do i = i1,i2 if (zeta(i) .LT. 0.0) then iq = iq + 1 iutb(iq) = i endif enddo if (iq .EQ. 0) go to 290 do i = 1,iq uzeta (i) = zeta(iutb(i)) ifz (i) = 1 uzetam(i) = 1.0e+50 ugzm (i) = 0.0e+00 enddo !------------------------------------------------------------------------ ! begin wegstein iteration for "zeta" at unstable points using ! hicks functions !------------------------------------------------------------------------ do icnt = 1,icntx do i = 1,iq if (ifz(i) .EQ. 0) go to 200 ugzzo = ALOG(zkmax(iutb(i))/ABS(zot(iutb(i)))) uzetao = ABS(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i) ux11 = 1. - 16.*uzeta(i) ux12 = 1. - 16.*uzetao y = SQRT(ux11) yo = SQRT(ux12) ufzo(i) = 1./yo ux13 = (1. + y)/(1. + yo) ux21 = ALOG(ux13) ufh(i) = (ugzzo - 2.*ux21)*rca x = SQRT(y) xo = SQRT(yo) xnum = (x**2 + 1.)*((x + 1.)**2) xden = (xo**2 + 1.)*((xo + 1.)**2) xtan = ATAN(x) - ATAN(xo) ux3 = ALOG(xnum/xden) uft(i) = (ugzzo - ux3 + 2.*xtan)*rca ! recompute scalers for ufm in terms of mom znot... zoc ugzzo = ALOG(zkmax(iutb(i))/ABS(zoc(iutb(i)))) uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i) ux11 = 1. - 16.*uzeta(i) ux12 = 1. - 16.*uzetao y = SQRT(ux11) yo = SQRT(ux12) ux13 = (1. + y)/(1. + yo) ux21 = ALOG(ux13) c ufzo(i) = 1./yo x = SQRT(y) xo = SQRT(yo) xnum = (x**2 + 1.)*((x + 1.)**2) xden = (xo**2 + 1.)*((xo + 1.)**2) xtan = ATAN(x) - ATAN(xo) ux3 = ALOG(xnum/xden) ufm(i) = (ugzzo - ux3 + 2.*xtan)*rca if (ngd .GE. nNEST-1) then !------------------------------------------------------------------------ ! obtain ten meter winds for diagnostic purposes !------------------------------------------------------------------------ ugz10 = ALOG(z10(iutb(i))/ABS(zoc(iutb(i)))) uzet1o = ABS(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i) uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i) ux11 = 1. - 16.*uzet1o ux12 = 1. - 16.*uzetao y = SQRT(ux11) y10 = SQRT(ux12) ux13 = (1. + y)/(1. + y10) ux21 = ALOG(ux13) x = SQRT(y) x10 = SQRT(y10) xnum = (x**2 + 1.)*((x + 1.)**2) xden = (x10**2 + 1.)*((x10 + 1.)**2) xtan = ATAN(x) - ATAN(x10) ux3 = ALOG(xnum/xden) uf10(i) = (ugz10 - ux3 + 2.*xtan)*rca endif ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i)) ux1 = (ugz - uzeta(i))/uzeta(i) ux2 = ABS(ux1) if (ux2 .GE. 5.0e-5) then uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i)) uzetam(i) = uzeta(i) if(uq .EQ. 1) then print*,'NCO: !!! ERROR !!! DIVIDE BY ZERO IN MFLUX2' & ,' (UNSTABLE CASE)' print*,'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i) endif uzeta (i) = (ugz - uzeta(i)*uq)/(1.0 - uq) ugzm (i) = ugz else ifz(i) = 0 endif 200 continue enddo enddo do i = 1,iq if (ifz(i) .GE. 1) go to 230 enddo go to 250 230 continue write(stdout,240) 240 format(2X, ' NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW ') call MPI_CLOSE(1,routine) !------------------------------------------------------------------------ ! gather unstable values !------------------------------------------------------------------------ 250 continue !------------------------------------------------------------------------ ! update "zo" for ocean points. zo cannot be updated within the ! wegsteins iteration as the scheme (for the near neutral case) ! can become unstable. !------------------------------------------------------------------------ do i = 1,iq uzo = zoc(iutb(i)) if (zoc(iutb(i)) .LT. 0.0) then wndm=wind(iutb(i))*0.01 ckg=0.0185*og uzo =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i)) cons_p000001 = .000001 cons_7 = 7. vis = 1.4E-1 ustar = sqrt( -uzo / zog) restar = -ustar * uzo / vis restar = max(restar,cons_p000001) ! Rat taken from Zeng, Zhao and Dickinson 1997 rat = 2.67 * restar ** .25 - 2.57 rat = min(rat ,cons_7) !constant rat=0.0 ! zot(iutb(i)) = uzo * exp(-rat) else zot(iutb(i)) = zoc(iutb(i)) endif !cbt zoc(iutb(i)) = uzo enddo !------------------------------------------------------------------------ ! obtain wind at ten meters for diagnostic purposes !------------------------------------------------------------------------ if (ngd .EQ. nNEST) then do i = 1,iq wind10(iutb(i),j) = uf10(i)*wind(iutb(i))/ufm(i) enddo endif do i = 1,iq xxfm(iutb(i)) = ufm(i) xxfm2(iutb(i)) = uft(i) xxfh(iutb(i)) = ufh(i) xxsh(iutb(i)) = ufzo(i) enddo 290 continue do i = i1,i2 ucom(i) = ukmax(i) vcom(i) = vkmax(i) if (windp(i) .EQ. 0.0) then windp(i) = 100.0 ucom (i) = 100.0/SQRT(2.0) vcom (i) = 100.0/SQRT(2.0) endif rho(i) = pss(i)/(rgas*(tsg(i) + enrca*(theta(i) - & tsg(i))*xxsh(i)/(xxfh(i) + enrca*xxsh(i)))) bq1(i) = wind(i)*rho(i)/(xxfm(i)*(xxfh(i) + enrca*xxsh(i))) enddo if (ntsflg .EQ. 0) go to 370 xks = 0.01 hcap = .5/2.39e-8 pith = SQRT(4.*ATAN(1.0)) alfus = alll/2.39e-8 teps = 0.1 !------------------------------------------------------------------------ ! pack land and sea ice points !------------------------------------------------------------------------ ip = 0 do i = i1,i2 if (land(i) .EQ. 1) then ip = ip + 1 indx (ip) = i slwa (ip) = - slwdc(i)/(2.39e-8*60.) tss (ip) = tstar(i) thetap (ip) = theta(i) rkmaxp (ip) = rkmax(i) aap (ip) = 5.673e-5 pssp (ip) = pss(i) ecofp (ip) = ecof(i) estsop (ip) = estso(i) rstsop (ip) = rstso(i) bq1p (ip) = bq1(i) bq1p (ip) = amax1(bq1p(ip),0.1e-3) delsrad(ip) = dt(1)*pith/(hcap*SQRT(3600.*24.*xks)) endif enddo !------------------------------------------------------------------------ ! initialize variables for first pass of iteration !------------------------------------------------------------------------ do i = 1,ip ifz (i) = 1 tsm (i) = tss(i) rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i))) if (nstep .EQ. -99 .AND. ngd .GT. 1 .OR. & nstep .EQ. -99 .AND. ngd .EQ. 1) then if (j .EQ. 1 .AND. i .EQ. 1) write(stdout,300) 300 format(2X, ' SURFACE EQUILIBRIUM CALCULATION ') foftm(i) = thetap(i) + 1./(cp*bq1p(i))*(slwa(i) - aap(i)* & tsm(i)**4 + ecofp(i)*alfus*bq1p(i)*rdiff(i)) else foftm(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsm(i)**4 - & cp*bq1p(i)*(tsm(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* & rdiff(i)) endif tsp(i) = foftm(i) enddo !------------------------------------------------------------------------ ! do iteration to determine "tstar" at new time level !------------------------------------------------------------------------ do icnt = 1,icntx do i = 1,ip if (ifz(i) .EQ. 0) go to 330 tab1 (i) = tsp(i) - 153.16 it (i) = IFIX(tab1(i)) tab2 (i) = tab1(i) - FLOAT(it(i)) t1 (i) = tab(it(i) + 1) t2 (i) = table(it(i) + 1) estsop(i) = t1(i) + tab2(i)*t2(i) psps2 = (pssp(i) - estsop(i)) if(psps2 .EQ. 0.0)then psps2 = .1 endif rstsop(i) = 0.622*estsop(i)/psps2 rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i))) if (nstep .EQ. -99 .AND. ngd .GT. 1 .OR. & nstep .EQ. -99 .AND. ngd .EQ. 1) then foft(i) = thetap(i) + (1./(cp*bq1p(i)))*(slwa(i) - aap(i)* & tsp(i)**4 + ecofp(i)*alfus*bq1p(i)*rdiff(i)) else foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - & cp*bq1p(i)*(tsp(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* & rdiff(i)) endif if (ngd .EQ. 1 .AND. j .EQ. 48 .AND. i .EQ. 19) then reflect = slwa(i) sigt4 = -aap(i)*tsp(i)**4 shfx = -cp*bq1p(i)*(tsp(i) - thetap(i)) alevp = ecofp(i)*alfus*bq1p(i)*rdiff(i) delten = delsrad(i) diffot = foft(i) - tss(i) endif frac(i) = ABS((foft(i) - tsp(i))/tsp(i)) !------------------------------------------------------------------------ ! check for convergence of all points use wegstein iteration !------------------------------------------------------------------------ if (frac(i) .GE. teps) then qf (i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i)) tsm (i) = tsp(i) tsp (i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i)) foftm(i) = foft(i) else ifz(i) = 0 endif 330 continue enddo enddo !------------------------------------------------------------------------ ! check for convergence of "t star" prediction !------------------------------------------------------------------------ do i = 1,ip if (ifz(i) .EQ. 1) then write(stdout, 340) tsp(i), i, j 340 format(2X, ' NON-CONVERGENCE OF T* PREDICTED (T*,I,J) = ', E14.8, & 2I5) write(stdout,345) indx(i), j, tstar(indx(i)), tsp(i), ip 345 format(2X, ' I, J, OLD T*, NEW T*, NPTS ', 2I5, 2E14.8, I5) write(stdout,350) reflect, sigt4, shfx, alevp, delten, diffot 350 format(2X, ' REFLECT, SIGT4, SHFX, ALEVP, DELTEN, DIFFOT ', & 6E14.8) call MPI_CLOSE(1,routine) endif enddo do i = 1,ip ii = indx(i) tstrc(ii) = tsp (i) enddo !------------------------------------------------------------------------ ! compute fluxes and momentum drag coef !------------------------------------------------------------------------ 370 continue do i = i1,i2 fxh(i) = bq1(i)*(theta(i) - tsg(i)) fxe(i) = ecof(i)*bq1(i)*(rkmax(i) - rstso(i)) if (fxe(i) .GT. 0.0) fxe(i) = 0.0 fxmx(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*ucom(i)/ & windp(i) fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ & windp(i) cdm(i) = 1./(xxfm(i)*xxfm(i)) cdm2(i) = 1./(xxfm2(i)*xxfm2(i)) enddo return end subroutine MFLUX2