! ---------------------------------------------------------------- subroutine mlo_main & (iyear, imonth, iday, & ilon, jlat, dz0, zlev0, lbot0, rr0, d10, d20, & coriolis0, qtot0, qsw0, taux0, tauy0, precip0, & evap0, qcor0, scor0, tclim0, sclim0, & fracice, h, t, s, laymix) ! ---------------------------------------------------------------- ! Fanglin Yang, July 2003 ! Purpose: 1-d mixed-layer ocean model, based on gaspar (jpo, 1988). ! note: ! 1. within the mixed layer, temperature and salinity are prognostics. ! mixed-layer depth is prognostics when entrainment occurs, and is ! diagnostics when mixed layer shoals. ! 2. in layers below the mixed layer, vertical diffusion, convective ! adjustment and solar radiation pentration are included. in ! addition, temperature and salinity are also relaxed to observed ! climatologies in a 10-year time scale. ! ! vertical structure of the model ! --------------------------- surface ! ! --------------------------- k=1, dz(1), zlev(1) ! ! --------------------------- k=2, dz(2), zlev(2) ! . ! . ! . ! --------------------------- ! ^^^ h, bottom of mixed layer ! --------------------------- k=laymix ! . ! . ! . ! --------------------------- k=km, dz(km), zlev(km) ! implicit none include 'mlo.h' ! input integer, intent(in) :: iyear !year integer, intent(in) :: imonth !month integer, intent(in) :: iday !day integer, intent(in) :: ilon !longitude index for 2-d application integer, intent(in) :: jlat !latitude index for 2-d application real, intent(in) :: dz0(km) !layer thickness, m real, intent(in) :: zlev0(km) !layer depth, m integer, intent(in) :: lbot0 !index for the bottom of the ocean real, intent(in) :: rr0 !ocean optical properties real, intent(in) :: d10 !ocean optical properties real, intent(in) :: d20 !ocean optical properties real, intent(in) :: coriolis0 !coriolis parameter f real, intent(in) :: qtot0 !surface downward total heat flux, w/m^2 real, intent(in) :: qsw0 !surface downward solar flux, w/m^2 real, intent(in) :: taux0 !surface zonal wind stress, n/m^2 real, intent(in) :: tauy0 !surface meridional wind stress, n/m^2 real, intent(in) :: precip0 !precipitation rate, kg/m^2/s real, intent(in) :: evap0 !evaporation rate, kg/m^2/s real, intent(in) :: qcor0 !heat flux correction term, w/m^2 real, intent(in) :: scor0 !salinity flux correction term, 1/1000 real, intent(in) :: tclim0(km) !temperature profile for relaxation, K real, intent(in) :: sclim0(km) !salinity profile for relaxation, K real, intent(in) :: fracice !ice fraction, 0-1 ! input & output real, intent(inout) :: h !mixed-layer depth real, intent(inout) :: t(km) !temperature profile, K real, intent(inout) :: s(km) !salinity profile, 1/1000 integer, intent(inout) :: laymix !layer in which the mixed-layer resides ! local real :: tmix, smix, taumag integer :: k !----------------------------------------------------------------------------- ! passing dummy variables that are defined in mlo.h do k=1,km dz(k)=dz0(k) zlev(k)=zlev0(k) tclim(k)=tclim0(k) sclim(k)=sclim0(k) enddo lbot =lbot0 rr =rr0 d1 =d10 d2 =d20 coriolis=coriolis0 qtot =qtot0 qsw =qsw0 precip=precip0 evap =evap0 qcor =qcor0 scor =scor0 taux =taux0 tauy =tauy0 ! ! insure a large enough u* which is obtained by assuming ! the wind speed being greater than ~ 1 m/sec equivalent ! to tau magnitude > 1.6e-3. ustar value under ice from ! fichefet and gaspar (1988 jpo pp 184) taumag = max(sqrt(taux**2 + tauy**2),0.001625) ustar = 0.006 * fracice + sqrt(taumag/rhoref)*(1-fracice) tmix = t(1) smix = s(1) do k=1, laymix-1 t(k) = tmix s(k) = smix end do ! ---------------------------------------------------------------- if(iwrite.eq.1) then write(99,*) "---------------- MLO_MAIN.F -------------------" write(99,*) "year,month,day,lon,lat",iyear, imonth, iday, ilon, jlat write(99,*) "dz(k),zlev(k),tclim(k),sclim(k),t(k),s(k)" do k=1,lbot write(99,'(i4,6f10.2)') k,dz(k),zlev(k),tclim(k),sclim(k),t(k),s(k) enddo write(99,*) "laymix, h, tmix, smix, fracice" write(99,'(i4,6f10.2)')laymix, h, tmix, smix, fracice write(99,*) "lbot, rr, d1, d2, coriolis" write(99,'(i4,3f10.2,e12.4)') lbot, rr, d1, d2, coriolis write(99,*)"qtot, qsw, taux, tauy, precip, evap" write(99,'(2f10.2, 4e12.4)')qtot, qsw, taux, tauy, precip, evap write(99,*)"ustar, qcor, scor" write(99,'(3f15.6)') ustar, qcor, scor endif ! ! ---------------------------------------------------------------- ! to represent the impact of currents and horizontal eddies ! in layers below the mixed layer, relaxtion is applied. call relax(tmix, smix, laymix, t, s) ! convective adjustment must be called before entrainment ! calculation to keep a stable stratification. this is because ! only when deltab>0, the assumption ap>0 <==> we>0 be true. call convect(h, tmix, smix, laymix, t, s) ! changes in temperature, salinity and mixed-layer depth due to entrainment call entrain(h, laymix, tmix, smix, t, s) ! changes in temperature and salinity due to ! surface heat and water flux, including flux corrections ! and penetration of solar radiation call surflux(h, tmix, smix, laymix, t, s) ! ! changes in temprature and salinity due to vertical diffusion call diffusion(h, laymix, km, lbot, iwrite, & visct, timestep, dz, zlev, t) call diffusion(h, laymix, km, lbot, iwrite, & viscs, timestep, dz, zlev, s) ! tmix=t(1) smix=s(1) if(iwrite.eq.1) write(99,*) "---------------- MLO_MAIN.F -------------------" do k=1,laymix if(t(k).gt.tmax .or. t(k).lt.tmin) then if(iwrite.eq.1) write(99,111) iyear, imonth, iday, ilon, jlat, k, t(k) write(6,111) iyear, imonth, iday, ilon, jlat, k, t(k) endif enddo do k=1,laymix if(s(k).gt.smax .or. s(k).lt.smin) then if(iwrite.eq.1) write(99,222) iyear, imonth, iday, ilon, jlat, k, s(k) write(6,222) iyear, imonth, iday, ilon, jlat, k, s(k) endif enddo 111 format("WARN t(k)", 6i6,f10.2) 222 format("WARN s(k)", 6i6,f10.2) end subroutine mlo_main ! ------------------------------------------------------------- subroutine apwe(h, tmix, smix, tbelow, sbelow, & alpha, beta, ap, flag, we) ! ------------------------------------------------------------- ! purpose: compute ap, and we if flag=1 ! following eq.(50) of gaspar, jpo, 1988. implicit none include 'mlo.h' real, parameter :: m1=0.45, m2=2.6, m3=1.9, m4=2.3, m5=0.6 real, parameter :: a1=0.6, a2=0.3 ! input real, intent(in) :: h, tmix, smix, tbelow, sbelow real, intent(in) :: alpha, beta integer, intent(in) :: flag ! output real, intent(out) :: ap real, intent(out) :: we ! local variables real :: bh, lamada,length,hl34,hlp48,sp,cp1,cp3,c4 real :: tmp1,tmp2,tmp3,tmp4, temp, deltab real :: rhomix, rhobelow ! ---------------------------------------------------------- if(iwrite.eq.1) then write(99,*) "-------------- APWE.F ------------------------- " write(99,*) "input: h tmix, smix, tbelow, sbelow, alpha, beta" write(99,'(5f12.6,2e12.4)') h, tmix, smix, tbelow, sbelow, alpha, beta endif call bh14(bh, h, tmix, smix, alpha, beta) !bouyancy lamada=abs(coriolis)/ustar !inverse of ekman lengthh scale length=bh/(ustar*ustar*ustar) !inverse of bulk monin-obukhov length temp=max(-30.,min(30.,h*length)) !provent from overflowing hl34=a1+a2*max(1.0,2.5*lamada*h)*exp(temp) !eq34 hlp48=a1+a2*exp(temp) !eq48 cp1=( (2.0-2.0*m5)*(hl34/hlp48) + m4 )/6.0 cp3=( m4*(m2+m3) - (hl34/hlp48)*(m2+m3-m5*m3) )/3.0 cp3=max(cp3,0.1) ap =cp3*ustar*ustar*ustar - cp1*h*bh ! compute we if(ap.gt.0.0 .and. flag.eq.1) then c4 =2.0*m4/(m1*m1) sp =(m2+m3)*ustar*ustar*ustar - 0.5*h*bh deltab=alpha*grav*(tmix-tbelow)-beta*grav*(smix-sbelow) !! add the parametrization of kim to the denominator to keep the !! mixed layer from over deepening deltab = deltab + 9.0*max(1.0e-4,ustar**2) tmp1=0.5*ap+cp1*sp tmp2=0.5*ap-cp1*sp tmp3=-tmp1+sqrt(tmp2*tmp2 + 2.0*c4*hl34*hl34*ap*sp) tmp4=c4*hl34*hl34 - cp1 we=tmp3/(tmp4*h*deltab) endif if(iwrite.eq.1) then write(99,*) "-------------- APWE.F ------------------------- " write(99,*) "output: ap, we, tmp1, tmp2, tmp3, tmp4, deltab" write(99,'(7e12.4)') ap, we, tmp1, tmp2, tmp3, tmp4, deltab endif end subroutine apwe !------------------------------------------------------ subroutine bh14(bh, h, tmix, smix, alpha, beta) !------------------------------------------------------ ! purpose: coumpte term B(h), eq.(14) of gaspar, jpo, 1988. implicit none include 'mlo.h' ! input real, intent(in) :: h !mixed-layer depth, meter real, intent(in) :: tmix !mixed-layer temperature, K real, intent(in) :: smix !mixed-layer salinity, 1/1000 real, intent(in) :: alpha !thermal expansion coefficient real, intent(in) :: beta !saline contration coefficient ! outpu real, intent(out) :: bh !bouyancy ! temporary real :: x1,x2,x3,x4,temp,sum,exph1,exph2 !------------------------------------------------------- if(iwrite.eq.1) then write(99,*) "------------- BH14.F ----------------" write(99,*) "input: h, tmix, smix, alpha, beta" write(99,'(3f10.2,2e12.4)') h, tmix, smix, alpha, beta endif ! determine solar penetration coefficient exph1=exp(max(-h/d1,-30.0)) exph2=exp(max(-h/d2,-30.0)) temp=alpha*grav/(rhoref*capa) x1=beta*grav*smix*(precip-evap)/rhoref x2=temp*qtot x3=temp*qsw*(rr*exph1+(1.0-rr)*exph2) ! vertical integral from the surface to depth (-h) sum=rr*d1*(1.0-exph1)+(1.0-rr)*d2*(1.0-exph2) x4=-2.0*temp*qsw*sum/h bh=x1+x2+x3+x4 if(iwrite.eq.1) then write(99,*) "output: bh, x1, x2, x3, x4" write(99,'(7e12.4)') bh, x1, x2, x3, x4 endif end subroutine bh14 ! --------------------------------------------------------- subroutine convect(h, tmix, smix, laymix, t, s) ! --------------------------------------------------------- ! purpose: determines if adjacent layers are convectively ! unstable, if they are, then mix them. mixed layer ! is included. implicit none include 'mlo.h' real, intent(inout) :: h real, intent(inout) :: tmix real, intent(inout) :: smix integer, intent(inout) :: laymix real, intent(inout) :: t(km) real, intent(inout) :: s(km) ! local real, allocatable :: dztmp(:),ttmp(:),stmp(:), rhotmp(:) real :: sumd, sumt, sums integer :: k, ktmp, ltop, llow logical :: flag_conv ! ---------------------------------------------- ! if(iwrite.eq.1) then write(99,*) "---------------- CONVECT.F --------------" write(99,*)" laymix, h, tmix, smix" write(99,'(i4, 3e12.4)') laymix, h, tmix, smix endif h=min(h,zlev(lbot)-0.5*dz(lbot)) flag_conv=.false. ! create new arrays including the mixed layer ktmp=lbot-laymix+2 allocate( dztmp(ktmp) ) allocate( ttmp(ktmp) ) allocate( stmp(ktmp) ) allocate( rhotmp(ktmp) ) dztmp(1)=h ttmp(1) =tmix stmp(1) =smix dztmp(2)=max(zlev(laymix)-h, 1.e-2) ttmp(2) =t(laymix) stmp(2) =s(laymix) if(ktmp.ge.3) then do k=3,ktmp dztmp(k)=dz(k+laymix-2) ttmp(k) =t(k+laymix-2) stmp(k) =s(k+laymix-2) enddo endif do k=1,ktmp call density(ttmp(k), max(stmp(k),0.0), rhotmp(k)) enddo ! convective adjustment do ltop=1,ktmp-1 do llow=ltop+1,ktmp if (rhotmp(ltop).gt.rhotmp(llow)) then if(ltop.eq.1) flag_conv=.true. sumd=0 sumt=0. sums=0. do k=ltop,llow sumt=sumt+ttmp(k)*dztmp(k)*rhotmp(k) sums=sums+stmp(k)*dztmp(k)*rhotmp(k) sumd=sumd+dztmp(k)*rhotmp(k) enddo do k=ltop,llow ttmp(k)=sumt/sumd stmp(k)=sums/sumd call density(ttmp(k), max(stmp(k),0.0), rhotmp(k)) enddo endif enddo enddo ! ! restore arrays tmix=ttmp(1) smix=stmp(1) do k=1,laymix-1 t(k)=ttmp(1) s(k)=stmp(1) enddo do k=laymix,lbot t(k)=ttmp(k-laymix+2) s(k)=stmp(k-laymix+2) enddo ! ! if convective adjustment does occur between the mixed-layer ! and layers below, find new laymix and h, which will be at ! the bottom edge of a layer where temperatue and salinity jumps occur. if(flag_conv) then laymix=1 do k=1,lbot-1 if(t(k).ne.tmix .or. s(k).ne.smix) goto 250 laymix=k+1 h=zlev(k)+1.0e-2 enddo 250 continue endif if(iwrite.eq.1) then write(99,*)"flag_conv=", flag_conv, " laymix, h, tmix, smix" write(99,'(i4, 3f10.2)') laymix, h, tmix, smix endif deallocate( dztmp ) deallocate( ttmp ) deallocate( stmp ) deallocate( rhotmp ) end subroutine convect ! ---------------------------------------- subroutine density(t, s, rho) ! ---------------------------------------- implicit none ! input real, intent(in) :: t !unit, K real, intent(in) :: s !unit, 1/1000 ! output real, intent(out) :: rho !unit, kg/m^3 ! local real :: tc ! compute density using the international equation ! of state of sea water 1980, (pond and pickard, ! introduction to dynamical oceanography, pp310). ! compression effects are not included rho = 0.0 tc = t - 273.15 ! effect of temperature on density (lines 1-3) ! effect of temperature and salinity on density (lines 4-8) rho = 1 999.842594 + 6.793952e-2 * tc 2 - 9.095290e-3 * tc**2 + 1.001685e-4 * tc**3 3 - 1.120083e-6 * tc**4 + 6.536332e-9 * tc**5 4 + 8.24493e-1 * s - 4.0899e-3 * tc * s 5 + 7.6438e-5 * tc**2 * s - 8.2467e-7 * tc**3 * s 6 + 5.3875e-9 * tc**4 * s - 5.72466e-3 * s**1.5 7 + 1.0227e-4 * tc * s**1.5 - 1.6546e-6 * tc**2 * s**1.5 8 + 4.8314e-4 * s**2 end subroutine density ! ---------------------------------------------------------------------------- subroutine diffusion(h, laymix, km, lbot, iwrite, & visc, dt, dz, zlev, u) ! ---------------------------------------------------------------------------- ! purpose: update temperature or salinity due to vertical diffusion implicit none ! input real, intent(in) :: h !mixed-layer depth integer, intent(in) :: laymix !layer in which mlo resides integer, intent(in) :: km !vertical layers integer, intent(in) :: lbot !ocean bottom layer integer, intent(in) :: iwrite !flag for diagnostic output real, intent(in) :: visc !viscocity real, intent(in) :: dt !timestep in seconds real, intent(in) :: dz(km) !layer thickness real, intent(in) :: zlev(km) !layer depth real, intent(inout) :: u(km) !temperature or salinity ! local variables integer :: k,j,kx,kmix real :: dx,r,aa,bb,cc,count,sum,unew(km) real, allocatable :: utmp(:), ztmp(:) real, allocatable :: y(:), f(:) real, allocatable :: gama(:), alpha(:), beta(:) ! ------------------------------------------------------------ ! ------------------------------------------------------------ ! METHOD DESCRIPTION ! the parabolic first-order partial differential equation ! which is of the typical form d(U)/d(t)-c*[d2(U)/d(x2)]=0, ! is sloved by classical implicit finite difference, ! [U(m,n+1)-U(m,n)] = r*[U(m+1,n+1)-2U(m,n+1)+U(m-1,n+1)] ! where r=c*dt/(dx*dx), m and n represent grids in space ! and time. m=1,2,...M-1, n=1,2,....N. ! ! For given ! initial condition U(m,n), m=0,1,2,....M, and ! boundary conditions U(0,n+1), U(M,n+1), ! to derive U(m,n+1), a linear tri-diagonal matix must be solved, ! | 1+2r -r |*| U(1,n+1) | = | U(1,n)+r*U(0,n+1) | ! | -r 1+2r -r |*| U(2,n+1) | = | U(2,n) | ! | -r 1+2r -r |*| U(3,n+1) | = | U(3,n) | ! | .............................. |*| ........ | = | ...... | ! | -r 1+2r |*| U(M-1,n+1)| = | U(M-1,n)+r*U(M,n+1)| ! ! Many standard routines can be used to solve this matrix. here I use the ! inverse-substitution method. The quantity in question is first ! projected to a uniform gird in vertical, and then restored back ! to oroginal non-uniform grid after diffusion calculation. ! ------------------------------------------------------------ ! ------------------------------------------------------------ ! ! search the smallest depth and set as the depth for uniform grid ! dx=dz(1) ! do k=2,lbot ! if(dz(k).lt.dx) dx=dz(k) ! enddo dx=5.0 !to save cpu, skip the search kx=int(zlev(lbot)/dx)-1 if(kx.lt.4) return ! allocate( utmp(0:kx+1) ) allocate( ztmp(0:kx+1) ) allocate( y(0:kx) ) allocate( f(0:kx) ) allocate( gama(0:kx) ) allocate(alpha(0:kx) ) allocate( beta(0:kx) ) r=visc*dt/(dx*dx) aa=-r bb=1.0+2.0*r cc=-r ! create profile in uniform vertical grid do k=0,kx+1 ztmp(k)=0.5*dx+k*dx !middle depth of each layer enddo kmix=int(h/dx)-1 do k=0,kmix !within the miex-layer utmp(k)=u(1) enddo do k=kmix+1,kx !below the miex-layer j=laymix-1 do while(ztmp(k).gt.zlev(j)) j=j+1 enddo utmp(k)=u(j) enddo ! set boundary conditions at k=0 and k=kx, right-hand-side matrix, ! determine coefficients for left-hand-side matrix for k=[1,kx-1] f(1)=utmp(1)+r*utmp(0) f(kx-1)=utmp(kx-1)+r*utmp(kx) do k=2,kx-2 f(k)=utmp(k) enddo alpha(1)=bb beta(1)=cc/alpha(1) do k=2,kx-1 gama(k)=aa alpha(k)=bb-gama(k)*beta(k-1) beta(k)=cc/alpha(k) enddo ! ! forward substitute y(1)=f(1)/bb do k=2,kx-1 y(k)=(f(k)-gama(k)*y(k-1))/alpha(k) enddo ! backward substitute utmp(kx-1)=y(kx-1) do k=kx-2,1,-1 utmp(k)=y(k)-beta(k)*utmp(k+1) enddo !----------------------------------------------- ! integrate and restore original non-uniform grid sum=0. count=0. do k=0,kmix !for mixed layer, derive mean u sum=sum+utmp(k) count=count+1.0 enddo do j=1,laymix-1 unew(j)=sum/count enddo k=kmix+1 sum=0. count=0. do j=laymix,lbot !below mixed layer do while(ztmp(k).lt.zlev(j)) sum=sum+utmp(k) count=count+1.0 k=k+1 enddo if(count.gt.0.0) then unew(j)=sum/count else unew(j)=u(j) !for case [zlev(laymix)]-h<0, mixed-layer deepens if(ap.gt.0.0 .and. h.lt.hmax ) then if(we.le.0.0) then ! note: in rare cases, when initialized from observations, ! such as in winter over Caspin Sea, the area ! is covered by ice, tmix=tice(271.36), smix is smaller ! than 10. then alpha<0 and beta>0. moreover, both temperature ! and salinity increase downward and hence rho increases downward. ! convective adjustment will not adjust the profiles. in this ! case, deltab<0 occurs and the condition ap>0 <==> we>0 is violated. ! entrainment is then skipped. if(iwrite.eq.1) then write(99,*) "inconsistence: ap>0, we<=0., skip" write(99,*) "ap, we, tmix, tbelow, smix, sbelow, alpha, beta" write(99,'(8e12.4)') ap, we, tmix, tbelow, smix, sbelow, alpha, beta endif goto 50 endif hnew=h+we*timestep if(hnew.gt.hmax) then hnew=hmax we=(hmax-h)/timestep endif ! determine new laymix after mixed layer deepens laynew = laymix do while (hnew .gt. zlev(laynew)) laynew = laynew+1 enddo ! determine new tbelow and sbelow using mass-weighted ! mean from hnew to h if(laynew.eq.laymix) then !within the same layer allmass= hnew-h allheat= t(laymix)*(hnew-h) allsalt= s(laymix)*(hnew-h) elseif(laynew.eq.laymix+1) then !cross one layer allmass= (zlev(laymix)-h) & + (hnew-zlev(laymix)) allheat= t(laymix)*(zlev(laymix)-h) & + t(laynew)*(hnew-zlev(laymix)) allsalt= s(laymix)*(zlev(laymix)-h) & + s(laynew)*(hnew-zlev(laymix)) else !corss more than one layer allmass= (zlev(laymix)-h) allheat= t(laymix)*(zlev(laymix)-h) allsalt= s(laymix)*(zlev(laymix)-h) do k=laymix+1,laynew-1 allmass=allmass + zlev(k) allheat=allheat + t(k)*zlev(k) allsalt=allsalt + s(k)*zlev(k) enddo allmass=allmass +(hnew-zlev(laynew-1)) allheat=allheat + t(laynew)*(hnew-zlev(laynew-1)) allsalt=allsalt + s(laynew)*(hnew-zlev(laynew-1)) endif tbelow=allheat/allmass sbelow=allsalt/allmass ! compute new tmix and smix tmixnew=tmix+timestep*we*(tbelow-tmix)/hnew smixnew=smix+timestep*we*(sbelow-smix)/hnew if(iwrite.eq.1) then write(99,*)"-------------- ENTRAIN.F -------------------" write(99,*)"mlo deepens, we=",we write(99,*)"laymix, laynew, h, hnew, tmix, tmixnew, smix, smixnew" write(99,'(2i4, 6f12.5)')laymix, laynew, h, hnew, tmix, tmixnew, smix, smixnew write(99,*)"tbelow, sbelow" write(99,'(6f12.5)')tbelow, sbelow endif 50 continue endif ! ---------------------------------------------------------- ! ---------------------------------------------------------- if (ap.le.0 .and. h .gt. dz(1)) then ! mixed-layer shoals. hnew satisfies ap(hnew)=0. use ! 2-point guessing method to solve ap=0. usually the function ! converges after less than 5 iterations for an accuracy ! of abs(ap)-->1.e-18 or abs[x(k+1)-x(k)]<1.e-2 . ! x(k+1)= x(k)-[x(k)-x(k-1)]*f[x(k)] / [f(x(k))-f(x(k-1))] ! note: ! 1. newton-raphson method can be used too, but needs to ! deduce d(ap)/d(h), which is nontrival either. ! 2. h is set no shallower than dz(1). ! 3. ap has a magnitude of about 1.e-5 to 1.0e-10 ! 4. sometimes the change in h is unreasonably large within ! one timestep. therefore, hnew-h is set no larger than ! 0.5*dz(laymix), such that h decreases in many timesteps. hmin=max(dz(1)+0.01,h-0.5*dz(laymix)) count=0 xkm1=h !first guess xk =h-0.5*dz(laymix) !second guess call apwe(xkm1, tmix, smix, tbelow, sbelow, & alpha, beta, apkm1, 0, we) call apwe(xk, tmix, smix, tbelow, sbelow, & alpha, beta, apk, 0, we) 100 xkp1=xk-(xk-xkm1)*apk/(apk-apkm1) xkp1=min(h,max(dz(1), xkp1)) !force dz(1)< hnew 10., 5., 5., 5., 5., 5., 5., 5., 5., 10., > 10., 10., 10., 10., 20., 20., 20., 20., 20., 40., > 40., 40., 60., 60., 60.,100., 100.,100.,100.,100./ ! all variables are local to the 1-D mixed-layer model mlo_main.F ! integer, parameter :: km =30 !model vertical layers real, parameter :: timestep=86400.0 !timestep, seconds real, parameter :: grav =9.81 !gravity, kg/m/s^2 real, parameter :: capa =3950.0 !heat capacity of sea water real, parameter :: rhoref = 1024.438 !sea water reference density, kg/m^3 real, parameter :: tmin=2.68E+02 !normal minimal temp real, parameter :: tmax=3.11E+02 !normal max temp real, parameter :: smin=1.0 !normal minimal salt real, parameter :: smax=50. !normal maximum salt real, parameter :: visct=1.e-5 !viscocity for temperature diffusion real, parameter :: viscs=1.e-5 !viscocity for salt diffusion integer, parameter :: iwrite = 0 !1--> output for dubug, 0--> no ! --constants integer :: lbot !level of ocean botom real :: dz(km) !depth, from top down, m real :: zlev(km) !sum(dz(k)), top down, m real :: rr !water optical coefficient real :: d1 !water optical coefficient real :: d2 !water optical coefficient common /mlocnst / dz, zlev, rr, d1, d2, lbot ! --input to mlo real :: precip !precipitation rate, kg/m^2/s real :: evap !evaporation rate, kg/m^2/s real :: qtot !surface total downward heat flux, w/m^2 real :: qsw !surface net downward solar flux, w/m^2 real :: taux !surface wind stress, n/m^2 real :: tauy !surface wind stress, n/m^2 real :: ustar !frictional velocity, m/s real :: coriolis !coriolis parameter f real :: qcor !heat flux correction, w/m^2 real :: scor !water flux correction, 1/1000 real :: tclim(km) !temperature profile for relaxation, K real :: sclim(km) !salinity profile for relaxation, K common /mloinput/ precip, evap, qtot, qsw, & taux, tauy, ustar, coriolis, & qcor, scor, tclim, sclim ! ---------------------------------------------- subroutine relax (tmix, smix, laymix, t, s) ! ---------------------------------------------- ! relax temperature and salinity below mixed layer to ! observed climatologies with a time scale of 10 years implicit none include 'mlo.h' real, parameter :: rate=3.170979198376e-9 !1/(10*365*86400) real, intent(inout) :: t(km) real, intent(inout) :: s(km) real, intent(inout) :: tmix real, intent(inout) :: smix integer, intent(in) :: laymix integer :: k do k=1, laymix-1 ! t(k) = t(k) + rate* (tclim(1) - t(k)) * timestep s(k) = s(k) + rate* (sclim(1) - s(k)) * timestep end do ! tmix = t(1) smix = s(1) do k=laymix, lbot t(k) = t(k) + rate* (tclim(k) - t(k)) * timestep s(k) = s(k) + rate* (sclim(k) - s(k)) * timestep end do end subroutine relax subroutine rhocoef(t, s, rhoref, alpha, beta) ! compute thermal expansion coefficient (alpha) ! and saline contraction coefficient (beta) using ! the international equation of state of sea water ! (1980). Ref: pond and pickard, introduction to ! dynamical oceanography, pp310. ! note: compression effects are not included implicit none real, intent(in) :: t, s, rhoref real, intent(out) :: alpha, beta real :: tc tc = t - 273.15 alpha = & 6.793952e-2 & - 2.0 * 9.095290e-3 * tc + 3.0 * 1.001685e-4 * tc**2 & - 4.0 * 1.120083e-6 * tc**3 + 5.0 * 6.536332e-9 * tc**4 & - 4.0899e-3 * s & + 2.0 * 7.6438e-5 * tc * s - 3.0 * 8.2467e-7 * tc**2 * s & + 4.0 * 5.3875e-9 * tc**3 * s & + 1.0227e-4 * s**1.5 - 2.0 * 1.6546e-6 * tc * s**1.5 ! alpha = -alpha/rhoref beta = & 8.24493e-1 - 4.0899e-3 * tc & + 7.6438e-5 * tc**2 - 8.2467e-7 * tc**3 & + 5.3875e-9 * tc**4 - 1.5 * 5.72466e-3 * s**.5 & + 1.5 * 1.0227e-4 * tc * s**.5 & - 1.5 * 1.6546e-6 * tc**2 * s**.5 & + 2.0 * 4.8314e-4 * s beta = beta / rhoref end subroutine rhocoef ! ------------------------------------------------------ subroutine surflux(h, tmix, smix, laymix, t, s) ! ------------------------------------------------------ ! purpose: update temperature and salinity due to ! surface heat and water flux, including flux corrections implicit none include 'mlo.h' ! input real, intent(in) :: h !mixed-layer depth, m integer, intent(in) :: laymix !layer in which mixed layer resides real, intent(inout) :: tmix !mixed-layer temperature, K real, intent(inout) :: smix !mixed-layer salinity, 1/1000 real, intent(inout) :: t(km) !temperature profile, K real, intent(inout) :: s(km) !salinity profile, 1/1000 ! local variables real :: qswpen, qswin, qswout, exph1, exph2, dist, dtk(km) real :: tnew(km), snew(km), tmixnew, smixnew integer :: k ! ---------------------------------------------------------- ! determine solar radiation that penetrates through the mixed layer exph1 =exp(max(-h/d1,-30.0)) exph2 =exp(max(-h/d2,-30.0)) qswpen=qsw*(rr*exph1+(1.0-rr)*exph2) ! mixed-layer temperature and salinity ! changes due to surface fluxes tmixnew=tmix + timestep*(qtot-qswpen+qcor)/(rhoref*capa*h) smixnew=smix + timestep*smix*(evap-precip+scor)/(rhoref*h) if(iwrite.eq.1) then write(99,*)"---------------- SURFLUX.F ----------------" write(99,*)"Changes by Surface FLuxes: h, tmix, Dtmix, smix, Dsmix" write(99,'(6f12.5)')h, tmix, tmixnew-tmix, smix, smixnew-smix endif tmix=tmixnew smix=smixnew do k=1,laymix-1 t(k)=tmix s(k)=smix enddo ! ------------------------------------------------------------ ! temperature changes in layers below the mixed-layer ! due to the penetration of solar radiation do k=1,km dtk(k)=0. enddo qswin=qswpen do 10 k=laymix,lbot exph1 =exp(max(-zlev(k)/d1,-30.0)) exph2 =exp(max(-zlev(k)/d2,-30.0)) qswout=qsw*(rr*exph1+(1.0-rr)*exph2) dist=dz(k) if(k.eq.laymix) dist=zlev(laymix)-h if(dist.le.1.e-2) goto 10 dtk(k)=timestep*(qswin-qswout)/(rhoref*capa*dist) t(k)=t(k) + dtk(k) qswin=qswout 10 continue if(iwrite.eq.1) then write(99,*)"---- t(k) change by solar penetration --------" write(99,*)"k, t(k), dtk(k)" do k=laymix, lbot write(99,'(i4, 2e13.5)')k, t(k)-dtk(k), dtk(k) enddo endif end subroutine surflux