module physics use const use comm use var use io implicit real*8 (a-h,o-z) !---------------------------------------------------------------------- ! namelist !---------------------------------------------------------------------- logical, parameter :: lmkcof_eos = .true. ! flag for making local_coeff logical, parameter :: eqstat_eos = .true. ! flag for the same eos as eq of state logical, parameter :: liceform = .true. ! flag for ice formatiom heat flux calc logical :: liceform_ini = .true. logical, parameter :: lactive_ice = .true. ! coupled to an active ice model real*8, parameter :: time_weight = 1.0 logical, parameter :: lcfc = .false. ! flag for cfc module logical :: lcfc_ini = .true. logical, parameter :: lkppmix = .true. ! flag for kpp vertical mixing scheme, parameters below are used logical :: lkpp_ini = .true. logical, parameter :: lrich = .true. ! compute Ri-dependent mixing logical, parameter :: ltidal_mixing = .true. ! tidal mixing logical, parameter :: lccsm_control_compatible = .false. logical, parameter :: lhoriz_varying_bckgrnd = .true. ! horizontal varying backgrd vdc,vvc logical, parameter :: ldbl_diff = .true. ! dbl diffusion logical, parameter :: lcheckekmo = .false. ! check Ekman, Monin-Obhukov depth limit logical, parameter :: lshort_wave = .true. ! computing short-wave forcing logical, parameter :: linertial = .false. ! inertial mixing parameterization logical, parameter :: llangmuir = .false. ! Langmuir parameterization logical, parameter :: smagmix = .false. ! flag for smagorinsky mixing scheme real*8, parameter :: hmix_tracer_type_gm = 1.0 !GM parameterization on !real*8, parameter :: hmix_tracer_type_gm = 0.0 !GM parameterization off logical, parameter :: ladjust_salinity_budget = .true. !---------------------------------------------------------------------- ! default coefficients and constants !---------------------------------------------------------------------- real*8, parameter :: convect_diff = 10000.d0, &! diffusivity to mimic convection convect_visc = 10000.d0 ! viscosity to mimic convection real*8, parameter :: epssfc = .1d0 ! non-dimensional extent of sfc layer !** parameters for Ri-dependent mixing integer, parameter :: num_v_smooth_Ri = 1 ! # of times to vertically smooth Ri real*8, parameter :: & Prandtl = 10.0, & rich_mix = 50.d0, &! coefficient for rich number term Riinfty = .8d0, &! Rich. no. limit for shear instability BVSQcon = .0d0 ! Brunt-Vaisala square cutoff(s**-2) !** parameters for dbl diffusion real*8, parameter :: & Rrho0 = 2.55d0, &! limit for double-diff density ratio dsfmax = 1.d0 ! max diffusivity for salt fingering !** parameters for bldepth real*8, parameter :: & vonkar = .4d0, &! von Karman constant cmonob = 1.d0, &! coefficient for Monin-Obukhov depth cekman = .7d0, &! coefficient for Ekman depth concv = 1.7d0 ! minimum allowed value for buoyancy frequency ratio at entrainment depth real*8, parameter, dimension(k1) :: & Ricr = .3d0 ! crit Rich no. for bndy layer depth ! real*8, parameter, dimension(i0,j0) :: bolus_sp = 0.d0 !scaled eddy-induced (bolus) speed used in inertial mixing !** parameters for velocity scale function (from Large et al.) real*8, parameter :: & zeta_m = -0.2d0, & zeta_s = -1.d0, & c_m = 8.38d0, & c_s = 98.96d0, & a_m = 1.26d0, & a_s = -28.86d0 !** parameters for subroutine blmix: mixing within boundary layer real*8, parameter :: & cstar = 10.d0 ! coeff for nonlocal transport !** parameters for tidal mixing real*8, parameter :: tidal_mix_max = 100.d0 !** grid info for mixing calculation real*8,dimension(0:k0) :: zgrid,hwide !** background vvc and vdc real*8, dimension(i0,j0,k1) :: bckgrnd_vvc, bckgrnd_vdc !** mixing within boundary layer real*8 :: & cg, & !coefficient for counter-gradient term Vtc !resolution and buoyancy independent part of the !turbulent velocity shear coefficient (for bulk Ri no) !** Langmuir parameterization real*8, dimension(i0,j0) :: fstokes !ratio of stokes velocity to ustar !** coriolis coeff real*8, dimension(j0) :: fcort !** tidal mixing real*8, dimension(i0,j0,k1) :: tidal_diff real*8, dimension(i0,j0,k1) :: tidal_coef contains ! ---------------------------------------------------------------------- subroutine windglo(u,v,days,dt,odz) dimension u(i0,j0),v(i0,j0),odz(*) !common/klimatglo/fnew,fold,new,nld !common/windsglo/taux(i2,j2,12),tauy(i2,j2,12) ! the velocity arrays u,v have a perimeter "ghost zone". ! thus, we set wind forcing only in the interior zones. ! taux,tauy are surface wind stress components. ! taux,tauy units are force per unit area (i.e., energy per unit volume). ! tmp=odz(1)/rho ! all units are cgs, so we use rho=1. if(myid.eq.0) write(*,*) 'days = ',days tmp=dt*odz(1) do j=2,j1 do i=2,i1 !u(i,j)=u(i,j)+tmp*(fold*taux(i-1,j-1,nld)+fnew*taux(i-1,j-1,new)) !v(i,j)=v(i,j)+tmp*(fold*tauy(i-1,j-1,nld)+fnew*tauy(i-1,j-1,new)) u(i,j)=u(i,j)+tmp*tauxx(i-1,j-1) v(i,j)=v(i,j)+tmp*tauyy(i-1,j-1) enddo enddo end subroutine ! ---------------------------------------------------------------------- subroutine cfcsurf(c,t,s,days,dt,odz) dimension c(i0,j0),t(i0,j0),s(i0,j0),odz(*) integer, parameter :: cfc_nt = 2 integer, parameter :: nyrbeg = 31 integer, parameter :: nyrend = 99 integer,parameter :: imt=i2 integer,parameter :: jmt=j2 real*8,save :: year(nyrend),p11(nyrend,cfc_nt),p12(nyrend,cfc_nt) real*8,save :: fice(i2,j2,12),xkw(i2,j2,12),patm(i2,j2,12) real*8,save :: ratio,day0,nday real*8,save :: p11_new(cfc_nt),p12_new(cfc_nt) real*8,save :: pnorth(cfc_nt),psouth(cfc_nt) real*8 :: rlon(imt,jmt),rlat(imt,jmt),CFCatm(imt,jmt,cfc_nt) real*8 :: Kw,CFCsat real*8 :: cfc_sc,cfc_alpha real*8, save :: CMAX,CMIN,STMAX,STMIN,SSMAX,SSMIN real*8, save :: KwMAX,scMAX,CFCMAX,alphaMAX,xkwMAX real*8, save :: KwMIN,scMIN,CFCMIN,alphaMIN,xkwMIN real*8 :: tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9 real*8 :: tmp10,tmp11,tmp12,tmp13,tmp14,tmp15,tmp16 INTEGER :: I,J,K,kn integer,save :: indx1,indx2,month !local if(lcfc_ini) then ! if(myid.eq.0)write(*,*)'days,dt=',days,dt if(myid.eq.0)write(*,*)'first call Read CFC gas and atm!' if(myid.eq.0)write(*,*)'i0,i1,i2=',i0,i1,i2 if(myid.eq.0)write(*,*)'j0,j1,j2=',j0,j1,j2 call read_gas(fice,xkw,patm) call read_cfcatm(cfc_nt,nyrbeg,nyrend,year,p11,p12) !Initial year in file indx1 = 90 !90 !31 !45 !1 !90 !we can use 29 =>spin up 2 yrs and run 37yrs indx2 = indx1+1 month=1 day0 = 0. !10800. if(myid.eq.0)write(*,*)'initial year,month:',1900+indx1,month STMAX = -9999.0 SSMAX = -9999.0 CMAX = -9999.0 KwMAX = -9999.0 scMAX = -9999.0 CFCMAX = -9999.0 alphaMAX = -9999.0 xkwMAX = -9999.0 STMIN = 9999.0 SSMIN = 9999.0 CMIN = 9999.0 KwMIN = 9999.0 scMIN = 9999.0 CFCMIN = 9999.0 alphaMIN = 9999.0 xkwMIN = 9999.0 do j=1,j2 do i=1,i2 rlon(i,j) = xdeg(i+1) rlat(i,j) = ydeg(j+1) enddo enddo lcfc_ini=.false. endif !first call nday = days-day0; if (mod(days,30.).eq.0.) then month = month +1 month = min(month,12) if(myid.eq.0)write(*,*)'month=',month end if if (mod(days,360.).eq.0.) then month = 1 day0 = days indx1 = indx2 indx2 = indx1+1 if(index2.ge.99) index2=99 !data only extend to 1999 if(myid.eq.0)write(*,*)'New year: month,year1,year2=',month,indx1,indx2 end if ratio = nday/360. !360.!360.0 p11_new(1) = p11(indx1,1)*(1-ratio)+p11(indx2,1)*ratio p11_new(2) = p11(indx1,2)*(1-ratio)+p11(indx2,2)*ratio p12_new(1) = p11(indx1,1)*(1-ratio)+p11(indx2,1)*ratio p12_new(2) = p11(indx1,2)*(1-ratio)+p11(indx2,2)*ratio pnorth(1)=p11_new(1) pnorth(2)=p12_new(1) psouth(1)=p11_new(2) psouth(2)=p12_new(2) call cfc_interp(cfc_nt,pnorth, psouth, rlat, imt, jmt, CFCatm) ! if(mod(days,1.).eq.0) then ! do j=1,J2 ! do i=1,I2 ! if(myid.eq.0)write(99,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1) ! if(myid.eq.1)write(98,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1) ! if(myid.eq.2)write(97,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1) ! if(myid.eq.3)write(96,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1) ! end do ! end do ! endif kn = 11 !11: CFC-11 12:CFC-12 ! if(myid.eq.0) write(*,*)'T,S=',T(70,60),S(70,60) ! if(myid.eq.0) call sc_cfc(T(70,60),kn,cfc_sc) ! if(myid.eq.0) call sol_cfc(T(70,60),S(70,60),kn,cfc_alpha) ! if(myid.eq.0)write(*,*)'alpha=',cfc_sc,cfc_alpha do j=1,j2 do i=1,i2 call sc_cfc(t(I+1,J+1),kn,cfc_sc) call sol_cfc(t(I+1,J+1),s(I+1,J+1),kn,cfc_alpha) CFCsat = cfc_alpha*CFCatm(I,J,1) !1=CFC-11 ! convert to pmol/kg for comparison with GLODAP data set CFCsat = CFCsat*1e12/1030. !pmol/kg Kw = (1-fice(I,J,month))*xkw(I,J,month)*(cfc_sc/660.0)**(-0.5) !m/s Kw = Kw*ODZ(1)*100. !the latter *100 because ODZ is in 1/cm rather than 1/m Kw = Kw*dt ! if(myid.eq.0) then ! if(i.eq.50.and.j.eq.50) then ! write(*,*)'check kw',myid,month,xkw(I+1,J+1,month) ! endif ! endif c(I+1,J+1) = c(I+1,J+1) + Kw*(CFCsat-c(I+1,J+1))*IN(I+1,J+1,1) STMAX=MAX(STMAX,T(I+1,J+1)) SSMAX=MAX(SSMAX,S(I+1,J+1)) CMAX=MAX(CMAX,C(I+1,J+1)) KwMAX=MAX(KwMAX,Kw) scMAX=MAX(scMAX,cfc_sc) CFCMAX=MAX(CFCMAX,CFCsat) alphaMAX=MAX(alphaMAX,cfc_alpha) xkwMAX=MAX(xkwMAX,xkw(I,J,month)) STMIN=MIN(STMIN,T(I+1,J+1)) SSMIN=MIN(SSMIN,S(I+1,J+1)) CMIN=MIN(CMIN,C(I+1,J+1)) KwMIN=MIN(KwMIN,Kw) scMIN=MIN(scMIN,cfc_sc) CFCMIN=MIN(CFCMIN,CFCsat) alphaMIN=MIN(alphaMIN,cfc_alpha) xkwMIN=MIN(xkwMIN,xkw(I,J,month)) end do end do tmp1=cmax tmp2=cmin tmp3=cfcmax tmp4=cfcmin tmp5=kwmax tmp6=kwmin tmp7=scmax tmp8=scmin tmp9=alphamax tmp10=alphamin tmp11=stmax tmp12=stmin tmp13=ssmax tmp14=ssmin tmp15=xkwmax tmp16=xkwmin CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp1,cmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp2,cmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp3,cfcmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp4,cfcmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp5,kwmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp6,kwmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp7,scmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp8,scmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp9,alphamax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp10,alphamin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp11,stmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp12,stmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp13,ssmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp14,ssmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp15,xkwmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp16,xkwmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) if(myid.eq.0) then if (mod(days,1.0).eq.0.) then write(*,*)'CFC check! DT=',DT write(*,*)days,month,nday,day0,ratio,indx1,indx2 write(*,*)stmax,ssmax,cmax,cfcmax,kwmax,scmax,alphamax write(*,*)stmin,ssmin,cmin,cfcmin,kwmin,scmin,alphamin write(*,*)xkwmax,xkwmin,odz(1),odz(1)*100.0*dt end if endif end subroutine subroutine sc_cfc(t,kn,cfc_sc) !--------------------------------------------------- ! CFC 11 and 12 schmidt number ! as a fonction of temperature. ! ! ref: Zheng et al (1998), JGR, vol 103,No C1 ! ! t: temperature (degree Celcius) ! kn: = 11 for CFC-11, 12 for CFC-12 ! ! J-C Dutay - LSCE !--------------------------------------------------- ! use cfc_glbl, only: rkind implicit none integer, intent(in) :: kn real*8 :: a1( 11: 12), a2( 11: 12), a3( 11: 12), a4( 11: 12) real*8,intent(in) :: t real*8,intent(out) :: cfc_sc ! ! coefficients with t in degre Celcius ! ------------------------------------ a1(11) = 3501.8 a2(11) = -210.31 a3(11) = 6.1851 a4(11) = -0.07513 ! a1(12) = 3845.4 a2(12) = -228.95 a3(12) = 6.1908 a4(12) = -0.067430 ! cfc_sc = a1(kn)+a2(kn)*t+a3(kn)*t*t+a4(kn)*t*t*t end subroutine subroutine sol_cfc(pt,ps,kn,cfc_alpha) implicit none real*8,intent(in) :: pt, ps real*8,intent(out) :: cfc_alpha real*8 :: ta,d real*8 :: a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12) real*8 :: b1 ( 11: 12), b2 ( 11: 12), b3 ( 11: 12) integer, intent(in) :: kn !! !! coefficient for solubility in mol/l/atm !! ---------------------------------------- ! ! for CFC 11 ! ---------- a1 ( 11) = -229.9261 a2 ( 11) = 319.6552 a3 ( 11) = 119.4471 a4 ( 11) = -1.39165 b1 ( 11) = -0.142382 b2 ( 11) = 0.091459 b3 ( 11) = -0.0157274 ! ! for CFC/12 ! ---------- a1 ( 12) = -218.0971 a2 ( 12) = 298.9702 a3 ( 12) = 113.8049 a4 ( 12) = -1.39165 b1 ( 12) = -0.143566 b2 ( 12) = 0.091015 b3 ( 12) = -0.0153924 ! ta = ( pt + 273.16)* 0.01 d = ( b3 ( kn)* ta + b2 ( kn))* ta + b1 ( kn) ! ! cfc_alpha = exp ( a1 ( kn) + a2 ( kn)/ ta + a3 ( kn)* log ( ta ) + & a4 (kn)* ta * ta + ps* d ) ! ! conversion from mol/(l * atm) to mol/(m^3 * atm) ! ------------------------------------------------ cfc_alpha = 1000. * cfc_alpha ! ! conversion from mol/(m^3 * atm) to mol/(m3 * pptv) ! -------------------------------------------------- cfc_alpha = 1.0e-12 * cfc_alpha end subroutine ! ---------------------------------------------------------------------- subroutine cfc_interp(cfc_nt,pnorth,psouth,rlat,imt,jmt,CFCatm) integer, parameter :: ijmax=200*200 integer, intent(in) :: cfc_nt, imt, jmt integer,save :: ientry integer :: i, j, ij, n real*8,save :: ys, yn real*8,save :: xintp(ijmax) real*8, intent(in) :: Pnorth(cfc_nt), Psouth(cfc_nt) real*8, intent(in) :: rlat(imt,jmt) real*8, intent(out) :: CFCatm(imt,jmt,cfc_nt) DATA ys, yn / -10., 10./ DATA ientry /0/ ientry = ientry + 1 ! Test to see if ijmax is large enough if (ijmax .lt. imt*jmt) then print *," cfc_interp: ERROR -> ijmax must be at least",imt*jmt stop endif ! IF Block to be executed ONLY on first call to this routine if (ientry .eq. 1) then do j = 1,jmt do i = 1,imt ij = (j-1)*imt + i if (rlat(i,j) .ge. yn) then xintp(ij) = 1. else if (rlat(i,j) .le. ys) then xintp(ij) = 0. else xintp(ij) = (rlat(i,j) - ys)/(yn-ys) endif end do end do endif ! Block to be executed every pass do n=1,cfc_nt do j=1,jmt do i=1,imt ij = (j-1)*imt + i CFCatm(i,j,n) = xintp(ij) * Pnorth(n) + & (1.0 - xintp(ij)) * Psouth(n) end do end do end do end subroutine ! ---------------------------------------------------------------------- subroutine read_cfcatm(cfc_nt,nyrbeg,nyrend,year,p11,p12) ! needs modification to dbl !local integer :: i, n integer, parameter :: iu = 10 integer,intent(in) :: nyrbeg,nyrend,cfc_nt real*8,intent(out) :: year(nyrend), p11(nyrend,cfc_nt), p12(nyrend,cfc_nt) character(len=80) :: filen year = 0. p11 = 0. p12 = 0. ! ! OPEN FILE ! --------- filen='prepglo/data/cfc1112_new.atm' open(UNIT=iu,FILE=filen,STATUS='old') ! ! Skip over 1st six descriptor lines ! ---------------------------------- do i=1,6 read(iu,99) end do ! ! READ FILE ! --------- do n = nyrbeg, nyrend READ(iu,*)year(n), p11(n,1), p12(n,1), p11(n,2), p12(n,2) year(n) = year(n) - 0.5 if(myid.eq.0) then write(*,*) 'read cfcatm file', year(n),p11(n,1), p12(n,1), p11(n,2), p12(n,2) endif end do close(iu) 99 FORMAT(1x) 100 FORMAT(f7.2, 4f8.2) end subroutine ! ---------------------------------------------------------------------- subroutine pp82glo ! ---------------------------------------------------------------------- ! vertical mixing coefficients ! approach: ! set background vertical diffusivities to dmz0 (o(0.1) cm-cm/sec) and ! add richardson number based vertical diffusivity plus numerical ! contribution according to vertical cell reynolds number. ! vertical cell re is o(100), except during winter cooling conditions ! and in wind-blown surface mixed layer, because internal waves, which ! dominate w below the sml in the model results during summer, do not ! mix t or s in nature. ! note: one cannot depend on diffusive closure by itself if one wants to ! model possible contra-diffusive effe! parameter(i01=i0+1) !common/windmxglo/vbk(k2),hbk(k2) !common/tauglo/add(i2,j2,k2) do k=1,k2 l=k+1 dzwl=orzmx/odzw(l) ! stability limit emax=.4/(dt*odzw(l)**2) do j=2,j1 do i=2,i1 ! pacanowski and philander (1981): evisc=a*(1/(1+b*ri))**n ! evisc=turbulent eddy viscosity; diffuse=turbulent diffusivity ! ri=max(0.,980.*(rho(i,j,l)-rho(i,j,k))*odzw(l)/(odzw(l)**2* ri=max(-0.02d0,980.*(rho(i,j,l)-rho(i,j,k))*odzw(l)/(odzw(l)**2*(0.001+(u2(i,j,l)-u2(i,j,k))**2+(v2(i,j,l)-v2(i,j,k))**2))) rrinv=1./(1.+5.*ri) tmp=rrinv**2 ! ri- and cell re- dependent vertical mixing ! numerically kosher because, for large ri, laminar flow dominates ! and time mean overshoot effe! are limited by time mean w being small temp=tmp*abs(w(i,j,l))*dzwl evisc=5.*tmp diffuse=evisc*rrinv ! add parameterizes extra momentum dissipation along steep slopes, ! high and low latitudes, and equatorial regions due to breaking of ! waves that were generated by big storm events that are missing when ! using climatological winds. add includes vbk(k) term. ev(i-1,j-1,k)=evisc+add(i-1,j-1,k)+temp hv(i-1,j-1,k)=diffuse+hbk(k)+temp ! apply stability limit and odzw factor tmp=odzw(l)*iw(i,j,l) ev(i-1,j-1,k)=tmp*min(emax,ev(i-1,j-1,k)) hv(i-1,j-1,k)=tmp*min(emax,hv(i-1,j-1,k)) enddo enddo enddo end subroutine ! ---------------------------------------------------------------------- subroutine pp82glonew ! ---------------------------------------------------------------------- ! vertical mixing coefficients ! approach: ! set background vertical diffusivities to dmz0 (o(0.1) cm-cm/sec) and ! add richardson number based vertical diffusivity plus numerical ! contribution according to vertical cell reynolds number. ! vertical cell re is o(100), except during winter cooling conditions ! and in wind-blown surface mixed layer, because internal waves, which ! dominate w below the sml in the model results during summer, do not ! mix t or s in nature. ! note: one cannot depend on diffusive closure by itself if one wants to ! model possible contra-diffusive effe! parameter(i01=i0+1) !common/windmxglo/vbk(k2),hbk(k2) !common/tauglo/add(i2,j2,k2) rate=-.05 jt=mylat*j2 do k=1,k2 l=k+1 dzwl=orzmx/odzw(l) ! deepen and augment mixing in high lats to emulate big storms tmpk=10.*exp(-.0001*z(2*k+1)) ! stability limit emax=.4/(dt*odzw(l)**2) elim=min(200.,emax) do j=2,j1 ! decrease background mixing except in high lats to emulate big ! storms within 20 grid intervals (~2 deg lat) of lat boundaries tmp=tmpk*(exp(rate*abs(2-j-jt)))+exp(rate*abs(j1t-j-jt)) vb=vbk(k)+tmp hb=hbk(k)+tmp do i=2,i1 ! pacanowski and philander (1981): evisc=a*(1/(1+b*ri))**n ! evisc=turbulent eddy viscosity; diffuse=turbulent diffusivity tmp=odzw(l) ! starting at day 0 of yr 16, decrease ri-based vertical mixing ! gap in izu ridge is also deepened to emulate unresolved deep kuroshio paths ! thus allowing it to pass over the ridgeline with less anomalous vortex ! squashing and thus taking a non-meander path when appropriate. oldfact=980.*(rho(i,j,l)-rho(i,j,k))*tmp/(tmp**2*(0.001+(u2(i,j,l)-u2(i,j,k))**2+(v2(i,j,l)-v2(i,j,k))**2)) ! big ri decreases vertical mixing ri=max(-0.19,5.*oldfact) rrinv=1./(1.+5.*ri) tmp=rrinv**2 ! ri- and cell re- dependent vertical mixing ! numerically kosher because, for large ri, laminar flow dominates ! and time mean overshoot effe! are limited by time mean w being small temp=.01*tmp*abs(w(i,j,l))*dzwl evisc=tmp diffuse=evisc*rrinv ! evisc, diffuse, hb and vb units are cm-cm/sec ev(i-1,j-1,k)=evisc+vb+temp hv(i-1,j-1,k)=diffuse+hb+0.1*temp ! apply stability limit and odzw factor tmp=odzw(l)*iw(i,j,l) ev(i-1,j-1,k)=tmp*min(elim,ev(i-1,j-1,k)) hv(i-1,j-1,k)=tmp*min(elim,hv(i-1,j-1,k)) enddo enddo enddo ! ====================================================================== ! add vertical mixing to emulate effect of equatorial vortices and ! penetrative convection in acc. ! ====================================================================== ! stacked vortices in 10-point equator band (hardwired to 1/4 deg res) ! strong velocity oscillations having ~20-day time scale are observed ! within a few deg of equator ! 33-point equatorial and acc bands ! ! do 770 i=1,i2 ! do 770 j=1,j2 ! do 770 k=1,k2 ! 770 scr(i+1,j+1,k)=ev(i,j,k) sc1=.01 sc2=.01 do k=1,k2 l=k+1 ! 100m scale height of added vertical mixing tmp=exp(-.0001*z(2*k+1)) do j=1,j2 tmp1=10./(1.+sc1*((0-ydeg(j))/(ydeg(j+1)-ydeg(j)))**2) addd=tmp1*odzw(l)*tmp do i=1,i2 ! insulated bottom, nonlinear drag invoked in fs ad=addd*iw(i+1,j+1,l) hv(i,j,k)=hv(i,j,k)+ad ev(i,j,k)=ev(i,j,k)+ad enddo enddo enddo end subroutine ! ---------------------------------------------------------------------- subroutine smagglo ! ---------------------------------------------------------------------- ! horizontal mixing coefficients ! approach:smagorinsky, j. general circulation experiments with the ! primitive equations, i. the basi! experiment. monthly weather rev. 1963, ! 91, 99-164. ! default value: 0.12 parameter(i01=i0+1) do k=1,k1 do j=2,j1 temp=ocs(j)*ody(j) do i=2,i1 scr(i,j,k)=1.2*dx(j)*dy(j)*in(i,j,k)*sqrt(((u(i,j,k)-u(i-1,j,k))*odx(j))**2+((csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)**2) enddo enddo enddo call mpi_sendrecv(scr(2,1,1),1,m_vlon,m_w,1,scr(i0,1,1),1,m_vlon,m_e,1,m_cart,istat,ierr) call mpi_sendrecv(scr(i1,1,1),1,m_vlon,m_e,1,scr(1,1,1),1,m_vlon,m_w,1,m_cart,istat,ierr) call mpi_sendrecv(scr(1,2,1),1,m_vlat,m_s,1,scr(1,j0,1),1,m_vlat,m_n,1,m_cart,istat,ierr) call mpi_sendrecv(scr(1,j1,1),1,m_vlat,m_n,1,scr(1,1,1),1,m_vlat,m_s,1,m_cart,istat,ierr) do k=1,k1 do j=1,j0 do i=1,i1 dmx(i,j,k)=dm0+.5*(scr(i,j,k)+scr(i+1,j,k)) enddo enddo do i=1,i0 dmy(i,1,k)=dm0+scr(i,2,k) dmy(i,j1,k)=dm0+scr(i,j2,k) do j=1,jf ! note: except in hi lats, dmy slightly violates ang. mom. conservation dmy(i,j,k)=dm0+.5*(scr(i,j,k)+scr(i,j+1,k)) enddo enddo enddo end subroutine ! ---------------------------------------------------------------------- subroutine salsurfglo(s) real*8, dimension(i0,j0), intent(inout) :: s real*8 dtmp ! Salinity= Salinity + [ tmp * sdot ] (sdot==evapo) ! tmp converts precipitation and evaporation rate (kg/m2/s) to psu ! hence tmp is DT/DZ /waterdensity * salinity ! tmp = dt*odz/1000 * psu ! from kg/m**2/s -> msu cm/s -> psu cm/s : -34.7*1.e-4*1.e3 dtmp=dt*odz(1)*1000.d0 !TSb dtmp=dtmp/2.0d0 DO j=1,j2 DO i=1,i2 ! S(I+1,J+1)=S(I+1,J+1)*(1-(DTMP*(EVAPO(i,j,1)+RAIN(I,J,1)))) + ! 1 DTMP*SALT(i,j,1) s(i+1,j+1)=s(i+1,j+1)+dtmp*STF(i+1,j+1,2) ! EVAPO<0 (negative up); PREC>0 (positive down); - (negative EVAPO+PREC) --> ! SALINITY INCREASE ! unit of SALT is kg(msu)/m^2/s enddo enddo end subroutine salsurfglo ! ---------------------------------------------------------------------- ! ---------------------------------------------------------------------- subroutine qsurfglo(t) real*8, dimension(i0,j0), intent(inout) :: t real*8 :: dtmp ! the temperature array t has a permiter "ghost zone". ! thus, we set atmospheri! heat exchange only in the interior zones. ! typical latitudinal range: -50 to 50 watts per square meter. ! 1 watt per square meter = 1 erg per second per equare cm. ! example: surface heating qdot = 50 watts per square meter !qdot=50. ! tmp is conversion factor from watts per square meter ! to deg ! per time step in top model layer. ! tmp=dt*odz(1)/rho/cp ! all units are cgs, so we use rho=1 !dog -> 1.026 ! cp for water is 2.5e4 ergs/gram/deg c. ! Cp for water is 3.996 !POP_CpSW = 3.996e7_POPd0 ! specific heat salt water !POP_CpSW = SHR_CONST_CPSW*10000.0_POPd0 ! erg/g/K dtmp=dt*odz(1)/(3.996d4*1.026) !TSb dtmp=dtmp/2.0d0 do j=1,j2 do i=1,i2 ! add because heat is going IN (positive downward) ! t(i+1,j+1)=t(i+1,j+1)+dtmp*(qdot(i,j,1)+qdot2(i,j,1)) ! t(i+1,j+1)=t(i+1,j+1)+dtmp*(qdot(i,j)+qdot2(i,j)) !TS t(i+1,j+1)=t(i+1,j+1)+dtmp*(qdot(i,j)+SHF_QSW(i+1,j+1)/hflux_factor) t(i+1,j+1)=t(i+1,j+1)+dtmp*qdot(i,j) enddo enddo end subroutine qsurfglo ! ---------------------------------------------------------------------- !subroutine sunglo(t,dt,odz,qprof,kb) ! integer*2 kb ! dimension t(i0,j0,k1),odz(k1),qprof(k1),kb(i0,j0) ! this is like subroutine qsurf except ! we use predetermined input array qprof to vertically distribute qdot. ! qprof(k) is the cumulative fraction of total heat absorbed to depth ! of bottom of layer k for a pure water column. ! example: surface radiation flux of 10 watts per square meter ! qdot=10. ! fract=qprof(1) ! do k=1,k1 ! dtemp=fract*qdot*dt/2.5e4 ! odz(k) is not the correct factor for layers deeper than the local bottom. ! thus, odz factor must be added inside loop 99. ! this quickly programmed logi! must be checked. ! do j=1,j2 ! do i=1,i2 ! kk=min(kb(i+1,j+1),k) ! t(i+1,j+1,kk)=t(i+1,j+1,kk)+dtemp*odz(kk) ! enddo ! enddo ! if (k.lt.k1) fract=qprof(k+1)-fract ! enddo !end subroutine !-------------------------------------------------------------------------- subroutine kppsrcglo(t,s) real*8 :: dtmp1, dtmp2 real*8, dimension(i0,j0,k1), intent(inout) :: t,s ! KPP_SRC unit: T (degree C/s) S (g/kg/S) ! SO ONLY NEED DT CONVERSION ! NEED TO CHECK THE UNIT OF STF/DZ*VDC again dtmp1 = dt dtmp2 = dt do k=1,k1 do j=1,j2 do i=1,i2 t(i+1,j+1,k)=t(i+1,j+1,k) + dtmp1*kpp_src(i+1,j+1,k,1)*IN(i+1,j+1,k) s(i+1,j+1,k)=s(i+1,j+1,k) + dtmp2*kpp_src(i+1,j+1,k,2)*IN(i+1,j+1,k) !add the KPP NON-LOCAL SOURCE enddo enddo enddo end subroutine kppsrcglo subroutine set_STF dtmp = (86400/daodt)*odz(1)*1000.d0 do j = 1,j2 do i = 1,i2 STF(i+1,j+1,1)=IN(i+1,j+1,1)*qdot(i,j)*hflux_factor ! W/m^2 to degC*cm/s STF(i+1,j+1,2)=IN(i+1,j+1,1)*((roff(i,j)+evapo(i,j)+prec(i,j)+& melt(i,j)+ioff(i,j))*(-34.7d0)*1.e-4+salt(i,j)*0.1d0)& !kg/m^2/s to msu*cm/s +(s_nudge(i+1,j+1)*in(i+1,j+1,1)*in(i+1,j+1,2)*& (fold*sclis(i+1,j+1,1,nld)+fnew*sclis(i+1,j+1,1,new)-s2(i+1,j+1,1)))/dtmp enddo enddo CALL MPI_EXCH_2D_R8(STF(:,:,1)) CALL MPI_EXCH_2D_R8(STF(:,:,2)) end subroutine set_STF subroutine kppglo(vdc) use forcing_coupled, only: QSW_COSZ_WGHT use blocks, only: nx_block, ny_block implicit none integer :: i,j,k,l real*8 :: emax,tmp real*8, dimension(i0,j0,k1) :: vvc real*8, dimension(i0,j0,0:k0,2), intent(inout) :: vdc ! !DESCRIPTION: ! This is the main driver routine which calculates EV,HV and KPP_SRC ! The KPP mixing scheme is outlined in Large et al, Reviews of Geophysics, 32, 363 (November 1994). ! KPP_SRC part is treated in ==surface effect== ! for testing, we use EP fluxes to reconstruct the missing STFs. ! swnet(:,:)=fold*qsw(:,:,nld)+fnew*qsw(:,:,new) ! qdot2=fold*qdot2tmp(:,:,nld)+fnew*qdot2tmp(:,:,new) ! qdot=fold*qdottmp(:,:,nld)+fnew*qdottmp(:,:,new) ! evapo=fold*evaptmp(:,:,nld)+fnew*evaptmp(:,:,new) ! prec=fold*prectmp(:,:,nld)+fnew*prectmp(:,:,new) ! use STF prec => deeper ! salt=fold*salttmp(:,:,nld)+fnew*salttmp(:,:,new) do j = 1,j2 do i = 1,i2 SHF_QSW(i+1,j+1)=IN(i+1,j+1,1)*qdot2(i,j)*hflux_factor ! W/m^2 to deg C*cm/s SMFT(i+1,j+1,1)=IN(i+1,j+1,1)*tauxx(i,j)!( fold*taux(i,j,nld) + fnew*taux(i,j,new) ) !dyne/cm^2 SMFT(i+1,j+1,2)=IN(i+1,j+1,1)*tauyy(i,j)!( fold*tauy(i,j,nld) + fnew*tauy(i,j,new) ) !dync/cm^2 enddo enddo SHF_QSW = QSW_COSZ_WGHT(2:nx_block-1,2:ny_block-1,1)*SHF_QSW CALL MPI_EXCH_2D_R8(SHF_QSW(:,:)) CALL MPI_EXCH_2D_R8(SMFT(:,:,1)) CALL MPI_EXCH_2D_R8(SMFT(:,:,2)) call qvolumeglo(t2,trans,0) !! flag = 0 ,calculate trans only call mpi_exch_3d_r8(u2) call mpi_exch_3d_r8(t2) call mpi_exch_3d_r8(s2) ! if (lopen==0) then ! call mpi_exch_3d_r8_vsymtry(v2) ! else call mpi_exch_3d_r8(v2) ! endif trcr(:,:,:,1)=t2 trcr(:,:,:,2)=s2*ppt_to_salt ! change to pop2 salt unit (g/g) call vmix_coeffs_kpp(IN,trans,kpp_src, vdc, vvc, trcr, u2,v2,rho,& stf, shf_qsw, convect_diff, convect_visc, & smft=smft) do k=1,k2 l=k+1 emax=.5d0/(dt*odzw(l)**2) !TS emax=MIN(1200.d0,emax) !emax=500.d0 do j=1,j2 do i=1,i2 ! if(k .ge. kbl(i+1,j+1))then ! emax=.5d0/(dt*odzw(l)**2) ! emax=MIN(6000.d0,emax) ! else ! emax=600.d0 ! endif dtmp=odzw(l)*iw(i+1,j+1,l) !! limit the value not to over 900 or .5/(DT*ODZW(L)**2) ev(i,j,k)=dtmp*MAX(MIN(emax,vvc(i+1,j+1,k)+add(i,j,k)),vbk(k)*5.d0) if ( hmix_tracer_type_gm < 0.5) then hv(i,j,k)=dtmp*MAX(MIN(emax,vdc(i+1,j+1,k,1)+add(i,j,k)),hbk(k)*5.d0) !if(myid==0) write(*,*) 'limited hv' endif end do end do end do end subroutine kppglo !-------------------------------------------------------------------------- subroutine vmix_coeffs_kpp(IN,trans,kpp_src, vdc, vvc, trcr, uuu,vvv,rhomix,& stf, shf_qsw, convect_diff, convect_visc, & smf,smft) ! !DESCRIPTION: ! This is the driver routine which calculates the vertical ! mixing coefficients for the KPP mixing scheme as outlined in ! Large, McWilliams and Doney, Reviews of Geophysics, 32, 363 ! (November 1994). The non-local mixing is also computed here, but ! is treated as a source term in baroclinic. implicit none ! !INPUT PARAMETERS: integer*2, dimension(i0,j0,k1), intent(in) :: IN real*8, dimension(i0,j0,k1,2), intent(inout) :: & trcr ! tracers at current time real*8, dimension(i0,j0,k1), intent(in) :: & uuu,vvv, & ! velocities at current time trans ! transmission (1 to 0) from surface (calculated in kppglo) real*8, dimension(i0,j0,k1), intent(in) :: & rhomix ! density at mix time real*8, dimension(i0,j0,2), intent(in) :: & stf ! surface forcing for all tracers real*8, dimension(i0,j0), intent(in) :: & shf_qsw ! short-wave forcing real*8, dimension(i0,j0,2), intent(in), optional :: & smf, &! surface momentum forcing at U points smft ! surface momentum forcing at T points ! *** either one or the other (not ! *** both) should be passed real*8, intent(in) :: & convect_diff, &! diffusivity to mimic convection convect_visc ! viscosity to mimic convection ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0,k1), intent(inout) :: & vvc ! viscosity for momentum diffusion real*8, dimension(i0,j0,0:k0,2),intent(inout) :: & vdc ! diffusivity for tracer diffusion real*8, dimension(i0,j0,k1,2), intent(out) :: & kpp_src !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & k, &! vertical level index i,j, &! horizontal loop indices n, &! tracer index mt2 ! index for separating temp from other trcrs integer, dimension(i0,j0) :: & kbl ! index of first lvl below hbl real*8, dimension(i0,j0) :: & ustar, &! surface friction velocity bfsfc, &! surface buoyancy forcing work1,work2,&! temporary storage fcon, &! convection temporary stable ! = 1 for stable forcing; = 0 for unstable forcing real*8, dimension(i0,j0,k1) :: & dbloc, &! buoyancy difference between adjacent levels dbsfc, &! buoyancy difference between level and surface ghat ! non-local mixing coefficient real*8, dimension(i0,j0,0:k0) :: & visc ! local temp for viscosity !----------------------------------------------------------------------- ! initialize !----------------------------------------------------------------------- if (lkpp_ini) then call init_vmix_kpp lkpp_ini=.false. endif !if(myid==73) write(*,*) 'trcrin2',trcr(1,1,1,2) !----------------------------------------------------------------------- ! compute buoyancy differences at each vertical level. ! derived variables: dbloc, dbsfc !----------------------------------------------------------------------- call buoydiff(dbloc, dbsfc, trcr) ! if(myid==73) write(*,*) 'trcrin3',trcr(1,1,1,2) !if (myid==73) write(*,*) 'dbloc=',dbloc(20,:,3) !if (myid==73) write(*,*) 'dbsfc=',dbsfc(20,:,3) !----------------------------------------------------------------------- ! compute mixing due to shear instability, internal waves and ! convection ! derived variables: vdc,visc !----------------------------------------------------------------------- call ri_iwmix(dbloc, visc, vdc, uuu, vvv) ! if(myid==73) write(*,*) 'trcrin4',trcr(1,1,1,2) ! if (myid==73) write(*,*) 'vdc=',vdc(20,:,4,1) ! if (myid==73) write(*,*) 'visc=',visc(20,:,4) !----------------------------------------------------------------------- ! compute double diffusion if desired !----------------------------------------------------------------------- if (ldbl_diff) call ddmix(vdc,trcr) ! if(myid==73) write(*,*) 'trcrin5',trcr(1,1,1,2) ! if (myid==73) write(*,*) 'ddmix=',vdc(20,:,4,1) !----------------------------------------------------------------------- ! compute boundary layer depth ! derived variables: kpp_hblt, ustar, bfsfc, stable, kbl !----------------------------------------------------------------------- if (present(smft)) then call bldepth (trans,dbloc, dbsfc, trcr, uuu, vvv, stf, shf_qsw, & kpp_hblt, ustar, bfsfc, stable, kbl, & smft=smft) else call bldepth (trans,dbloc, dbsfc, trcr, uuu, vvv, stf, shf_qsw, & kpp_hblt, ustar, bfsfc, stable, kbl, & smf=smf) endif ! if (myid==73) write(*,*) 'mld=',kpp_hblt(20,:) ! if (myid==73) write(*,*) 'uuu=',uuu(20,:,4) ! if (myid==73) write(*,*) 'stf11=',stf(20,:,1) ! if (myid==73) write(*,*) 'stf22=',stf(20,:,2) !----------------------------------------------------------------------- ! compute boundary layer diffusivities ! derived variables: visc,vdc,ghat !----------------------------------------------------------------------- call blmix(visc, vdc, kpp_hblt, ustar, bfsfc, stable, & kbl, ghat) !----------------------------------------------------------------------- ! ! consider interior convection: ! ! compute function of Brunt-Vaisala squared for convection. ! ! use either a smooth ! ! WORK1 = N**2, FCON is function of N**2 ! FCON = 0 for N**2 > 0 ! FCON = [1-(1-WORK1/BVSQcon)**2]**3 for BVSQcon < N**2 < 0 ! FCON = 1 for N**2 < BVSQcon ! ! or a step function. The smooth function has been used with ! BVSQcon = -0.2e-4_dbl_kind. ! ! after convection, average viscous coeffs to U-grid and reset sea ! floor values ! !----------------------------------------------------------------------- do k=1,k1-1 work1 = dbloc(:,:,k)/(zgrid(k) - zgrid(k+1)) if (BVSQcon /= 0.d0) then work2 = min(1.d0-(max(work1,BVSQcon))/BVSQcon, 1.d0) fcon = (1.d0 - work2*work2)**3 else where (work1 > 0.d0) fcon = 0.d0 elsewhere fcon = 1.d0 end where endif !if (myid==73) write(*,*) 'fcon',FCON(20,:),convect_visc,convect_diff !*** add convection and reset sea floor values to zero do j=1,j0 do i=1,i0 if ( k >= kbl(i,j) ) then visc(i,j,k) = visc(i,j,k) + convect_visc * fcon(i,j) vdc(i,j,k,1) = vdc(i,j,k,1) + convect_diff * fcon(i,j) vdc(i,j,k,2) = vdc(i,j,k,2) + convect_diff * fcon(i,j) endif if ( k >= kb(i,j) ) then visc(i,j,k ) = 0.d0 vdc (i,j,k,1) = 0.d0 vdc (i,j,k,2) = 0.d0 endif end do end do vvc(:,:,k) = merge(visc(:,:,k), 0.d0, ( k < kb(:,:) ) ) enddo VDC(:,:,k1,:) = 0. VVC(:,:,k1) = 0. !----------------------------------------------------------------------- ! add ghatp term from previous computation to right-hand-side ! source term on current row !----------------------------------------------------------------------- do n=1,2 mt2=min(n,2) KPP_SRC(:,:,1,n) = STF(:,:,n)*odz(1) & *(-VDC(:,:,1,mt2)*GHAT(:,:,1)) do k=2,k1 KPP_SRC(:,:,k,n) = STF(:,:,n)*odz(k) & *( VDC(:,:,k-1,mt2)*GHAT(:,:,k-1) & -VDC(:,:,k ,mt2)*GHAT(:,:,k )) end do end do KPP_SRC(:,:,:,2) = 1000.d0*KPP_SRC(:,:,:,2) !from msu to psu !if (myid==73) write(*,*) 'kppsrc1',kpp_src(20,:,1,1) ! if (myid==73) write(*,*) 'kppsrc2',kpp_src(20,:,1,2) !----------------------------------------------------------------------- ! ! compute diagnostic mixed layer depth (cm) using a max buoyancy ! gradient criterion. Use USTAR and BFSFC as temps. ! !----------------------------------------------------------------------- ustar = 0.d0 where (kb(:,:) == 1) hmxl(:,:) = z(2)-z(1) elsewhere hmxl(:,:) = 0.d0 endwhere do k=2,k1 where (k <= kb(:,:)) ustar = max(dbsfc(:,:,k)/(z(2*k)-z(1)),ustar) hmxl(:,:) = z(2*k)-z(1) endwhere enddo visc(:,:,1) = 0.d0 do k=2,k1 where (ustar > 0.d0 ) visc(:,:,k) = (dbsfc(:,:,k)-dbsfc(:,:,k-1))/ & (z(2*k) - z(2*k-2)) end where where ( visc(:,:,k) >= ustar .and. & (visc(:,:,k)-visc(:,:,k-1)) /= 0.d0 .and. & ustar > 0.d0 ) ! avoid divide by zero bfsfc = (visc(:,:,k) - ustar)/ & (visc(:,:,k)-visc(:,:,k-1)) hmxl(:,:) = -0.5d0*(zgrid(k ) + zgrid(k-1))*(1.0d0-bfsfc) & -0.5d0*(zgrid(k-1) + zgrid(k-2))*bfsfc ustar(:,:) = 0.d0 endwhere enddo end subroutine vmix_coeffs_kpp !-------------------------------------------------------------------------- subroutine qvolumeglo(t,transmitted,flag) ! tmp is conversion factor from watts per square meter ! to deg ! per time step in top model layer. ! tmp=dt*odz(1)/rho/cp ! all units are cgs, so we use rho=1 !dog -> 1.026 ! cp for water is 2.5e4 ergs/gram/deg c. ! Cp for water is 3.996 implicit none real*8, dimension(i0,j0,k1), intent(inout) :: t,transmitted integer, intent(in) :: flag !-------------------------------------------------------------------------- ! local variables !-------------------------------------------------------------------------- real*8 :: myZ,R,xi1,xi2,total_trans,dtmp real*8 :: qdot2sum,qdot2ori,qdotori integer :: i,j,k R=.58d0 xi1=.35d0 xi2=23.d0 do j=1,j0 do i=1,i0 total_trans=0.d0 do k=1,k1 myZ=-z(2*k)/100.d0 transmitted(i,j,k)=(R*exp(myZ/xi1)+(1-R)*exp(myZ/xi2))* & in(i,j,k)*odz(k) total_trans=total_trans+transmitted(i,j,k) enddo do k=1,k1 if (total_trans.gt.1.d-6) then transmitted(i,j,k)=transmitted(i,j,k)/total_trans else transmitted(i,j,k)=0.d0 endif enddo enddo enddo if (flag==0) goto 8888 ! absorption*QDOT is added to T !TS qdot2ori=SUM(qdot2(:,:)*in(2:i1,2:j1,1)) !TS dtmp=qdot2ori !TS CALL MPI_ALLREDUCE(DTMP,QDOT2ORI,1,MPI_REAL8,MPI_SUM,M_CART,IERR) !TS qdotori=SUM(qdot(:,:)*in(2:i1,2:j1,1)) !TS dtmp=qdotori !TS CALL MPI_ALLREDUCE(DTMP,QDOTORI,1,MPI_REAL8,MPI_SUM,M_CART,IERR) qdot2sum=0.d0 do k=1,K1 dtmp=dt*odz(k)/(3.996d4*1.026) !TSb dtmp=dtmp/2.0d0 do j=1,j2 do i=1,i2 t2(i+1,j+1,k)=t2(i+1,j+1,k) + & !TS dtmp*qdot2(i,j)*transmitted(i+1,j+1,k)*IN(i+1,j+1,k) !add because heat is going IN !TS qdot2sum=qdot2sum+qdot2(i,j)*transmitted(i+1,j+1,k)*IN(i+1,j+1,k) dtmp*SHF_QSW(i+1,j+1)/hflux_factor*transmitted(i+1,j+1,k)*IN(i+1,j+1,k) !add because heat is going IN qdot2sum=qdot2sum+SHF_QSW(i+1,j+1)/hflux_factor*transmitted(i+1,j+1,k)*IN(i+1,j+1,k) enddo enddo enddo !TS dtmp=qdot2sum !TS CALL MPI_ALLREDUCE(DTMP,QDOT2SUM,1,MPI_REAL8,MPI_SUM,M_CART,IERR) !TS if(myid==0)WRITE(*,*)'TOTALQ:',QDOTORI,QDOT2ORI,QDOT2SUM 8888 continue end subroutine qvolumeglo !-------------------------------------------------------------------------- !-------------------------------------------------------------- subroutine buoydiff(dbloc,dbsfc,trcr) ! !DESCRIPTION: ! This routine calculates the buoyancy differences at model levels. ! ! !REVISION HISTORY: ! same as module ! implicit none ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1,2), intent(inout) :: & trcr ! tracers at current time ! !OUTPUT PARAMETERS: real*8 , dimension(i0,j0,k1), intent(out) :: & dbloc, &! buoyancy difference between adjacent levels dbsfc ! buoyancy difference between level and surface !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & k, &! vertical level index i,j, &! horizontal indices kprev, klvl, ktmp ! indices for 2-level TEMPK array integer, dimension(i0,j0) :: kmt real*8 , dimension(i0,j0) :: & rho1, &! density of sfc t,s displaced to k rhokm, &! density of t(k-1),s(k-1) displaced to k rhok, &! density at level k tempsfc ! adjusted temperature at surface real*8 , dimension(i0,j0,2) :: & tempk ! temp adjusted for freeze at levels k,k-1 real*8, dimension(i0,j0) :: test,test2,outs !----------------------------------------------------------------------- ! calculate density and buoyancy differences at surface !----------------------------------------------------------------------- tempsfc = merge(-2.d0,trcr(:,:,1,1),trcr(:,:,1,1) < -2.d0) klvl = 2 kprev = 1 tempk(:,:,kprev) = tempsfc dbsfc(:,:,1) = 0.d0 !test = 24.d0 !test2 = 34.d0/1000.d0 !call state(4,4,test(:,:),test2(:,:),RHOFULL=outs(:,:) ) !if (myid==4) write(*,*) 'test',outs(5,5) !----------------------------------------------------------------------- ! calculate DBLOC and DBSFC for all other levels !-----------------------------------------------------------------------\ do k = 2,k1 tempk(:,:,klvl) = merge(-2.0d0,trcr(:,:,k,1),trcr(:,:,k,1) < -2.0d0) ! calculate the density from Temperature and Salinity using an equation of state ! derived from McDougall, Wright, Jackett and Feistel (hereafter MWJF, 2001 ) call state(k, k, tempsfc, trcr(:,:,1,2), & rhofull=rho1) call state(k, k, tempk(:,:,kprev), trcr(:,:,k-1,2), & rhofull=rhokm) call state(k, k, tempk(:,:,klvl), trcr(:,:,k ,2), & rhofull=rhok) !if (k==4 .and. myid==73) write(*,*) 'tempk',tempk(20,:,klvl) !if (k==4 .and. myid==73) write(*,*) 'trcr',trcr(20,:,k,2) do j=1,j0 do i=1,i0 if (rhok(i,j) /= 0.d0) then dbsfc(i,j,k) = grav*(1.0d0 - rho1 (i,j)/rhok(i,j)) dbloc(i,j,k-1) = grav*(1.0d0 - rhokm(i,j)/rhok(i,j)) else dbsfc(i,j,k) = 0.d0 dbloc(i,j,k-1) = 0.d0 endif if ( k-1 >= kb(i,j) ) dbloc(i,j,k-1) = 0.d0 enddo enddo ktmp = klvl klvl = kprev kprev = ktmp !if (k==4 .and. myid==73) write(*,*) 'dbcalc=',dbsfc(20,:,k) enddo dbloc(:,:,k1) = 0.d0 end subroutine buoydiff !-------------------------------------------------------------- !-------------------------------------------------------------- subroutine state(k, kk, tempk, saltk, & rhoout, rhofull, drhodt, drhods) ! !DESCRIPTION: ! Returns the density of water at level k from equation of state ! $\rho = \rho(d,\theta,S)$ where $d$ is depth, $\theta$ is ! potential temperature, and $S$ is salinity. the density can be ! returned as a perturbation (RHOOUT) or as the full density ! (RHOFULL). Note that only the polynomial EOS choice will return ! a perturbation density; in other cases the full density is returned ! regardless of which argument is requested. ! ! This routine also computes derivatives of density with respect ! to temperature and salinity at level k from equation of state ! if requested (ie the optional arguments are present). ! ! If $k = kk$ are equal the density for level k is returned. ! If $k \neq kk$ the density returned is that for a parcel ! adiabatically displaced from level k to level kk. ! implicit none ! !INPUT PARAMETERS: integer, intent(in) :: & k, &! depth level index kk ! level to which water is adiabatically ! displaced real*8, dimension(i0,j0), intent(in) :: & tempk, &! temperature at level k saltk ! salinity at level k ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0), optional, intent(out) :: & rhoout, &! perturbation density of water rhofull, &! full density of water drhodt, &! derivative of density with respect to temperature drhods ! derivative of density with respect to salinity !EOP !BOC !----------------------------------------------------------------------- ! local variables: !----------------------------------------------------------------------- integer :: dep real*8, dimension(i0,j0) :: & TQ,SQ, &! adjusted T,S sqr,denomk, &! work arrays work1,work2,work3,work4 real*8 :: ppz ! temporary pressure scalars real*8 :: & mwjfnums0t0, mwjfnums0t1, mwjfnums0t2, mwjfnums0t3, & mwjfnums1t0, mwjfnums1t1, mwjfnums2t0, & mwjfdens0t0, mwjfdens0t1, mwjfdens0t2, mwjfdens0t3, mwjfdens0t4, & mwjfdens1t0, mwjfdens1t1, mwjfdens1t3, & mwjfdensqt0, mwjfdensqt2 !----------------------------------------------------------------------- ! valid ranges and pressure as function of depth !----------------------------------------------------------------------- real*8, dimension(k1) :: & tmin, tmax, &! valid temperature range for level k smin, smax, &! valid salinity range for level k pressz ! ref pressure (bars) at each level !----------------------------------------------------------------------- ! MWJF EOS coefficients !----------------------------------------------------------------------- !*** these constants will be used to construct the numerator !*** factor unit change (kg/m^3 -> g/cm^3) into numerator terms real*8, parameter :: p001 = 0.001_r8 real*8, parameter :: & mwjfnp0s0t0 = 9.99843699e+2_r8 * p001, & mwjfnp0s0t1 = 7.35212840e+0_r8 * p001, & mwjfnp0s0t2 = -5.45928211e-2_r8 * p001, & mwjfnp0s0t3 = 3.98476704e-4_r8 * p001, & mwjfnp0s1t0 = 2.96938239e+0_r8 * p001, & mwjfnp0s1t1 = -7.23268813e-3_r8 * p001, & mwjfnp0s2t0 = 2.12382341e-3_r8 * p001, & mwjfnp1s0t0 = 1.04004591e-2_r8 * p001, & mwjfnp1s0t2 = 1.03970529e-7_r8 * p001, & mwjfnp1s1t0 = 5.18761880e-6_r8 * p001, & mwjfnp2s0t0 = -3.24041825e-8_r8 * p001, & mwjfnp2s0t2 = -1.23869360e-11_r8* p001 !*** these constants will be used to construct the denominator real*8, parameter :: & mwjfdp0s0t0 = 1.0e+0_r8, & mwjfdp0s0t1 = 7.28606739e-3_r8, & mwjfdp0s0t2 = -4.60835542e-5_r8, & mwjfdp0s0t3 = 3.68390573e-7_r8, & mwjfdp0s0t4 = 1.80809186e-10_r8, & mwjfdp0s1t0 = 2.14691708e-3_r8, & mwjfdp0s1t1 = -9.27062484e-6_r8, & mwjfdp0s1t3 = -1.78343643e-10_r8, & mwjfdp0sqt0 = 4.76534122e-6_r8, & mwjfdp0sqt2 = 1.63410736e-9_r8, & mwjfdp1s0t0 = 5.30848875e-6_r8, & mwjfdp2s0t3 = -3.03175128e-16_r8, & mwjfdp3s0t1 = -1.27934137e-17_r8 !*** prevent problems with garbage on land points or ghost cells TQ = min(tempk, 999.0d0) TQ = max(TQ, -2.0d0) SQ = min(saltk, 0.999d0) SQ = max(SQ, 0.d0) !------------------------------------------------------------ do dep=1,k1 pressz(dep) = pressure(-zgrid(dep)*.01_r8) ! cm to m enddo !if (myid==73) write(*,*) 'pressz',pressz(:) !if (myid==73) write(*,*) 'zt=',-zgrid(1:k1) ppz = 10.d0*pressz(kk) SQ = 1000.d0*SQ !#ifdef CCSMCOUPLED ! call shr_vmath_sqrt(SQ, SQR, nx_block*ny_block) !#else SQR = sqrt(SQ) !#endif !*** !*** first calculate numerator of MWJF density [P_1(S,T,p)] !*** mwjfnums0t0 = mwjfnp0s0t0 + ppz*(mwjfnp1s0t0 + ppz*mwjfnp2s0t0) mwjfnums0t1 = mwjfnp0s0t1 mwjfnums0t2 = mwjfnp0s0t2 + ppz*(mwjfnp1s0t2 + ppz*mwjfnp2s0t2) mwjfnums0t3 = mwjfnp0s0t3 mwjfnums1t0 = mwjfnp0s1t0 + ppz*mwjfnp1s1t0 mwjfnums1t1 = mwjfnp0s1t1 mwjfnums2t0 = mwjfnp0s2t0 WORK1 = mwjfnums0t0 + TQ * (mwjfnums0t1 + TQ * (mwjfnums0t2 + & mwjfnums0t3 * TQ)) + SQ * (mwjfnums1t0 + & mwjfnums1t1 * TQ + mwjfnums2t0 * SQ) !*** !*** now calculate denominator of MWJF density [P_2(S,T,p)] !*** mwjfdens0t0 = mwjfdp0s0t0 + ppz*mwjfdp1s0t0 mwjfdens0t1 = mwjfdp0s0t1 + ppz**3 * mwjfdp3s0t1 mwjfdens0t2 = mwjfdp0s0t2 mwjfdens0t3 = mwjfdp0s0t3 + ppz**2 * mwjfdp2s0t3 mwjfdens0t4 = mwjfdp0s0t4 mwjfdens1t0 = mwjfdp0s1t0 mwjfdens1t1 = mwjfdp0s1t1 mwjfdens1t3 = mwjfdp0s1t3 mwjfdensqt0 = mwjfdp0sqt0 mwjfdensqt2 = mwjfdp0sqt2 WORK2 = mwjfdens0t0 + TQ * (mwjfdens0t1 + TQ * (mwjfdens0t2 + & TQ * (mwjfdens0t3 + mwjfdens0t4 * TQ))) + & SQ * (mwjfdens1t0 + TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3)+ & SQR * (mwjfdensqt0 + TQ*TQ*mwjfdensqt2)) DENOMK = 1.d0/WORK2 if (present(RHOOUT)) then RHOOUT = WORK1*DENOMK endif if (present(RHOFULL)) then RHOFULL = WORK1*DENOMK endif if (present(DRHODT)) then WORK3 = &! dP_1/dT mwjfnums0t1 + TQ * (2.d0*mwjfnums0t2 + & 3.d0*mwjfnums0t3 * TQ) + mwjfnums1t1 * SQ WORK4 = &! dP_2/dT mwjfdens0t1 + SQ * mwjfdens1t1 + & TQ * (2.d0*(mwjfdens0t2 + SQ*SQR*mwjfdensqt2) + & TQ * (3.d0*(mwjfdens0t3 + SQ * mwjfdens1t3) + & TQ * 4.d0*mwjfdens0t4)) DRHODT = (WORK3 - WORK1*DENOMK*WORK4)*DENOMK endif if (present(DRHODS)) then WORK3 = &! dP_1/dS mwjfnums1t0 + mwjfnums1t1 * TQ + 2.d0*mwjfnums2t0 * SQ WORK4 = mwjfdens1t0 + &! dP_2/dS TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3) + & 1.5d0*SQR*(mwjfdensqt0 + TQ*TQ*mwjfdensqt2) DRHODS = (WORK3 - WORK1*DENOMK*WORK4)*DENOMK * 1000.d0 !unit converison endif end subroutine state function pressure(depth) ! !INPUT PARAMETERS: real*8, intent(in) :: depth ! depth in meters ! !OUTPUT PARAMETERS: real*8 :: pressure ! pressure [bars] pressure = 0.059808d0*(exp(-0.025d0*depth) - 1.0d0) & + 0.100766d0*depth + 2.28405e-7*depth**2 end function pressure subroutine ri_iwmix(dbloc, visc, vdc, uuu, vvv) ! !DESCRIPTION: ! Computes viscosity and diffusivity coefficients for the interior ! ocean due to shear instability (richardson number dependent), ! internal wave activity, and to static instability (Ri < 0). implicit none ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1), intent(in) :: & uuu ! U velocities at current time real*8, dimension(i0,j0,k1), intent(in) :: & vvv, &! V velocities at current time dbloc ! buoyancy difference between adjacent levels ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0,0:k0,2), intent(inout) :: & vdc ! diffusivity for tracer diffusion ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0,0:k0), intent(inout) :: & visc ! viscosity !EOP !BOC !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & k, &! index for vertical levels i,j, &! horizontal loop indices n ! vertical smoothing index real*8, dimension(i0,j0) :: & vshear, &! (local velocity shear)^2 ri_loc, &! local Richardson number fri, &! function of Ri for shear work1 real*8, dimension(i0,j0) :: kvmix,kvmix_m real*8, dimension(i0,j0,0:k0) :: work0 !----------------------------------------------------------------------- ! compute mixing at each level !----------------------------------------------------------------------- kvmix = 0.d0 kvmix_m = 0.d0 work0(:,:,0) = 0.d0 do k = 1,k1 !----------------------------------------------------------------------- ! compute velocity shear squared and average to T points: ! VSHEAR = (UUU(k)-UUU(k+1))**2+(VVV(k)-VVV(k+1))**2 ! Use FRI here as a temporary. !----------------------------------------------------------------------- if (k < k1) then fri = (uuu(:,:,k)-uuu(:,:,k+1))**2 + & (vvv(:,:,k)-vvv(:,:,k+1))**2 vshear=fri else vshear = 0.d0 endif !----------------------------------------------------------------------- ! compute local richardson number ! use visc array as temporary Ri storage to be smoothed !----------------------------------------------------------------------- ri_loc = dbloc(:,:,k)*(zgrid(k)-zgrid(k+1))/(vshear + eps) work0(:,:,k) = merge( ri_loc, work0(:,:,k-1), k <= int(kb(:,:)) ) enddo !----------------------------------------------------------------------- ! vertically smooth Ri num_v_smooth_Ri times with 1-2-1 weighting ! result again stored temporarily in WORK0 and use RI_LOC and FRI ! as temps !----------------------------------------------------------------------- do n = 1,num_v_smooth_Ri fri = 0.25d0 * work0(:,:,1) work0(:,:,k1+1) = work0(:,:,k1) do k=1,k1 do j=1,j0 do i=1,i0 ri_loc(i,j) = work0(i,j,k) if ( int(kb(i,j)) >= 3 ) then !modified from KMT, land=0 work0(i,j,k) = fri(i,j) + 0.5d0*ri_loc(i,j) & + 0.25d0*work0(i,j,k+1) endif fri(i,j) = 0.25d0*ri_loc(i,j) enddo enddo enddo enddo !----------------------------------------------------------------------- ! now that we have a smoothed Ri field, finish computing coeffs ! at each level !----------------------------------------------------------------------- if ( ltidal_mixing ) tidal_diff(:,:,:) = 0.d0 ! write(*,*) 'kilmer_ri3' do k = 1,k1 !----------------------------------------------------------------------- ! ! if Ri-number mixing requested, ! evaluate function of Ri for shear instability: ! for 0 < Ri < Riinfty, function = (1 - (Ri/Riinfty)**2)**3 ! for Ri > Riinfty, function = 0 ! for Ri < 0 , function = 1 ! compute contribution due to shear instability ! WORK0 holds smoothed Ri at k, but replaced by real VISC ! ! otherwise only use iw ! convection is added later ! !----------------------------------------------------------------------- if ( ltidal_mixing ) then !----------------------------------------------------------------------- ! ! consider the internal wave mixing first. rich_mix is used as the ! upper limit for internal wave mixing coefficient. bckgrnd_vvc ! was already multiplied by Prandtl. ! ! NOTE: no partial_bottom_cell implementation at this time ! !----------------------------------------------------------------------- work1 = dbloc(:,:,k)/(zgrid(k) - zgrid(k+1)) where (work1 > 0.d0) tidal_diff(:,:,k) = tidal_coef(:,:,k)/work1 endwhere ! Notes: ! (1) this step breaks backwards compatibility ! (2) check for k>2 was added to if statement to ! avoid out of bounds access if (.not. lccsm_control_compatible) then ! this step breaks backwards compatibility if ( k > 2 ) then where ( ( k ==kb(:,:)-1 .or. k == kb(:,:)-2 ) ) tidal_diff(:,:,k) = max( tidal_diff(:,:,k), & tidal_diff(:,:,k-1)) endwhere endif endif work1 = Prandtl*min(bckgrnd_vvc(:,:,k)/Prandtl & + tidal_diff(:,:,k), tidal_mix_max) if ( k < k1 ) then kvmix_m(:,:) = work1(:,:) endif if ( k < k1 ) then vdc(:,:,k,2) = min(bckgrnd_vdc(:,:,k) + tidal_diff(:,:,k), & tidal_mix_max) kvmix(:,:) = vdc(:,:,k,2) endif if (lrich) then fri = min((max(work0(:,:,k),0.d0))/Riinfty, 1.d0) visc(:,:,k) = work1 + rich_mix*(1.d0 - fri*fri)**3 if ( k < k1 ) then vdc(:,:,k,2) = vdc(:,:,k,2) + rich_mix*(1.d0 - fri*fri)**3 vdc(:,:,k,1) = vdc(:,:,k,2) endif else visc(:,:,k) = work1 if ( k < k1 ) then vdc(:,:,k,1) = vdc(:,:,k,2) endif endif else ! .not. ltidal_mixing ! write(*,*) 'kilmer_km1',k if ( k < k1 ) then kvmix(:,:) = bckgrnd_vdc(:,:,k) kvmix_m(:,:) = bckgrnd_vvc(:,:,k) endif if (lrich) then fri = min((max(work0(:,:,k),0.d0))/Riinfty, 1.0d0) visc(:,:,k ) = bckgrnd_vvc(:,:,k) + & rich_mix*(1.0d0 - fri*fri)**3 if ( k < k1 ) then vdc(:,:,k,2) = bckgrnd_vdc(:,:,k) + & rich_mix*(1.0d0 - fri*fri)**3 vdc(:,:,k,1) = vdc(:,:,k,2) endif else visc(:,:,k ) = bckgrnd_vvc(:,:,k) if ( k < k1 ) then vdc(:,:,k,2) = bckgrnd_vdc(:,:,k) vdc(:,:,k,1) = vdc(:,:,k,2) endif endif endif ! ltidal_mixing !----------------------------------------------------------------------- ! set seafloor values to zero !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 if ( k >= int(kb(i,j)) ) then visc(i,j,k ) = 0.d0 vdc (i,j,k,1) = 0.d0 vdc (i,j,k,2) = 0.d0 endif enddo enddo enddo !----------------------------------------------------------------------- ! fill extra coefficients for blmix !----------------------------------------------------------------------- visc(:,:,0 ) = 0.d0 vdc (:,:,0,:) = 0.d0 visc(:,:,k1+1 ) = 0.d0 vdc (:,:,k1+1,:) = 0.d0 end subroutine ri_iwmix subroutine init_vmix_kpp implicit none !----------------------------------------------------------------------- ! local coefficients !----------------------------------------------------------------------- integer :: i,j,k real*8 :: bckgrnd_vdc_psis ! PSI diffusivity in northern hemisphere real*8 :: bckgrnd_vdc_psin ! PSI diffusivity in southern hemisphere real*8, dimension(k1) :: zw ! vert dist from sfc to bottom of layer real*8, parameter :: bckgrnd_vdc1 = 0.16d0 ! background diffusivity (Ledwell) real*8, parameter :: bckgrnd_vdc2 = 0.0 ! variation in diffusivity real*8, parameter :: bckgrnd_vdc_eq = 0.01 ! equatorial diffusivity (Gregg) real*8, parameter :: bckgrnd_vdc_psim = 0.13 ! Max. PSI induced diffusivity (MacKinnon) real*8, parameter :: bckgrnd_vdc_ban = 1.0 ! Banda Sea diffusivity (Gordon) real*8, parameter :: bckgrnd_vdc_dpth = 1000.0e02 !2500 ! depth at which diff equals vdc1 real*8, parameter :: bckgrnd_vdc_linv = 4.5e-05 ! inverse length for transition region ! tidal mixing real*8, dimension(i0,j0) :: & vertical_func, & work1,HT real*8 :: coef,rho_fw real*8, parameter :: & local_mixing_fraction = 0.33, & mixing_efficiency = 0.2, & vertical_decay_scale = 500.0e02, & tidal_mix_max = 100.0 if (lhoriz_varying_bckgrnd) then k = 1 do j=1,j0 do i=1,i0 bckgrnd_vdc_psis= bckgrnd_vdc_psim*exp(-(0.4d0*(ydeg(j)+28.9d0))**2.) bckgrnd_vdc_psin= bckgrnd_vdc_psim*exp(-(0.4d0*(ydeg(j)-28.9d0))**2.) bckgrnd_vdc(i,j,k)=bckgrnd_vdc_eq+bckgrnd_vdc_psin+bckgrnd_vdc_psis if ( ydeg(j) .lt. -10.d0 ) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc(i,j,k) + bckgrnd_vdc1 elseif ( ydeg(j) .le. 10.d0 ) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc(i,j,k) + & bckgrnd_vdc1 * (ydeg(j)/10.d0)**2. else bckgrnd_vdc(i,j,k)=bckgrnd_vdc(i,j,k) + bckgrnd_vdc1 endif !---------------- ! North Banda Sea !---------------- if ( (ydeg(j) .lt. -1.0d0) .and. (ydeg(j) .gt. -4.0d0) .and. & (xdeg(i) .gt. 103.0d0) .and. (xdeg(i) .lt. 134.0d0)) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc_ban endif !----------------- ! Middle Banda Sea !----------------- if ( (ydeg(j) .le. -4.0d0) .and. (ydeg(j) .gt. -7.0d0) .and. & (xdeg(i) .gt. 106.0d0) .and. (xdeg(i) .lt. 140.d0)) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc_ban endif !---------------- ! South Banda Sea !---------------- if ( (ydeg(j) .le. -7.0d0) .and. (ydeg(j) .gt. -8.3d0) .and. & (xdeg(i) .gt. 111.0d0) .and. (xdeg(i) .lt. 142.0d0)) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc_ban endif bckgrnd_vvc(i,j,k) = Prandtl*bckgrnd_vdc(i,j,k) end do ! i end do ! j do k=2,k1 bckgrnd_vdc(:,:,k) = bckgrnd_vdc(:,:,1) bckgrnd_vvc(:,:,k) = bckgrnd_vvc(:,:,1) enddo else ! .not. lhoriz_varying_bckgrnd do k=1,k1 zw(k) = z(2*k+1) bckgrnd_vdc(:,:,k) = bckgrnd_vdc1 + bckgrnd_vdc2* & atan(bckgrnd_vdc_linv* & (zw(k)-bckgrnd_vdc_dpth)) bckgrnd_vvc(:,:,k) = Prandtl*bckgrnd_vdc(:,:,k) if (bckgrnd_vdc2 /= 0.d0 .and. myid == 0) then write (*,*)'bckgrnd_vdc=', bckgrnd_vdc(1,1,k) endif end do endif ! lhoriz_varying_bckgrnd ! initialize grid info hwide(0) = eps zgrid(0) = eps ! zgrid(1) = -0.5d0/odz(1) zgrid(1) = -z(2) do k =1,k1-1 !zgrid(k+1) = - (-zgrid(k) + 1.d0/odzw(k+1)) !zgrid(k+1) = -z(2*(k+1)) zgrid(k+1) = zgrid(k)-(z(2*k+1)-z(2*k-1) + z(2*k+3)-z(2*k+1))/2 enddo do k=1,k1 !hwide(k) = 1.d0/odz(k) hwide(k) = z(2*k+1)-z(2*k-1) enddo zgrid(k0) = -z(k0+k1) hwide(k0) = eps ! initialize cg and Vtc !if (myid==0) then !write(*,*) 'zgrid=',zgrid !write(*,*) 'hwide=',hwide !endif Vtc = sqrt(0.2d0/c_s/epssfc)/vonkar**2 cg = cstar*vonkar*(c_s*vonkar*epssfc)**0.33d0 ! stokes ratio do j = 1,j0 fstokes(:,j) = 11.0d0 - max(5.0d0*cos(3.0d0* (ydeg(j)/180.d0*pi) ) ,0.d0 ) enddo ! coriolis coeff do j=1,j0 fcort(j) = 2.d0*7.292123625e-5_r8*sin(ydeg(j)/180.d0*pi) enddo if (ltidal_mixing) then tidal_diff = 0.d0 tidal_coef = 0.d0 rho_fw = 1.0d0 HT = 0.0 do k = 1,k1 do j = 1,j0 do i = 1,i0 if (k == int(kb(i,j))) HT(i,j) = z(2*k+1) enddo enddo enddo !! tidal_energy now move to ./src/driver.F90 for model initialization. !! (see. set windmix) !----------------------------------------------------------------------- ! convert TIDAL_ENERGY_FLUX from W/m^2 to gr/s^3. !----------------------------------------------------------------------- ! tidal_energy_flux(:,:) = 1000.d0 * tidal_energy_flux(:,:) !----------------------------------------------------------------------- ! compute the time independent part of the tidal mixing coefficients !----------------------------------------------------------------------- work1=0.d0 do k=1,k1 where ( k < int(kb(:,:)) ) work1(:,:) = work1(:,:) + & exp(-(HT(:,:) - z(2*k+1))/vertical_decay_scale) / odzw(k+1) !dzw(k) endwhere end do ! k tidal_coef(:,:,:) = 0.d0 coef = local_mixing_fraction * mixing_efficiency / rho_fw do k=1,k1 where ( k <= int(kb(:,:)) ) vertical_func(:,:) = & exp(-(HT(:,:) - z(2*k+1)) & /vertical_decay_scale) / work1(:,:) tidal_coef(:,:,k) = coef * & tidal_energy_flux(:,:) * vertical_func(:,:) elsewhere tidal_coef(:,:,k) = 0.d0 endwhere enddo ! k call mpi_exch_3d_r8(tidal_coef) end if ! ltidal_mixing end subroutine init_vmix_kpp subroutine ddmix(VDC, TRCR) ! !DESCRIPTION: ! $R_\rho$ dependent interior flux parameterization. ! Add double-diffusion diffusivities to Ri-mix values at blending ! interface and below. implicit none ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1,2), intent(in) :: & TRCR ! tracers at current time ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0,0:k0,2),intent(inout) :: & VDC ! diffusivity for tracer diffusion !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: k,kup,knxt real*8, dimension(i0,j0) :: & ALPHADT, &! alpha*DT across interfaces BETADS, &! beta *DS across interfaces RRHO, &! dd density ratio DIFFDD, &! dd diffusivity scale prtl_tmp ! prandtl number real*8, dimension(i0,j0,2) :: & TALPHA, &! temperature expansion coeff SBETA ! salinity expansion coeff !----------------------------------------------------------------------- ! compute alpha*DT and beta*DS at interfaces. use RRHO and ! PRANDTL for temporary storage for call to state !----------------------------------------------------------------------- kup = 1 knxt = 2 prtl_tmp = merge(-2.0d0,TRCR(:,:,1,1),TRCR(:,:,1,1) < -2.0d0) call state(1, 1, prtl_tmp, TRCR(:,:,1,2),RHOFULL=RRHO, & DRHODT=TALPHA(:,:,kup), DRHODS=SBETA(:,:,kup)) do k=1,k1 if ( k < k1 ) then prtl_tmp = merge(-2.0d0,TRCR(:,:,k+1,1),TRCR(:,:,k+1,1) < -2.0d0) call state(k+1, k+1, prtl_tmp, TRCR(:,:,k+1,2), & RHOFULL=RRHO, DRHODT=TALPHA(:,:,knxt), & DRHODS= SBETA(:,:,knxt)) ALPHADT = -0.5d0*(TALPHA(:,:,kup) + TALPHA(:,:,knxt)) & *(TRCR(:,:,k,1) - TRCR(:,:,k+1,1)) BETADS = 0.5d0*( SBETA(:,:,kup) + SBETA(:,:,knxt)) & *(TRCR(:,:,k,2) - TRCR(:,:,k+1,2)) kup = knxt knxt = 3 - kup else ALPHADT = 0.d0 BETADS = 0.d0 endif !----------------------------------------------------------------------- ! salt fingering case !----------------------------------------------------------------------- where ( ALPHADT > BETADS .and. BETADS > 0.d0 ) RRHO = MIN(ALPHADT/BETADS, Rrho0) DIFFDD = dsfmax*(1.0d0-(RRHO-1.0d0)/(Rrho0-1.0d0))**3 VDC(:,:,k,1) = VDC(:,:,k,1) + 0.7d0*DIFFDD VDC(:,:,k,2) = VDC(:,:,k,2) + DIFFDD endwhere !----------------------------------------------------------------------- ! diffusive convection !----------------------------------------------------------------------- where ( ALPHADT < 0.d0 .and. BETADS < 0.d0 .and. ALPHADT > BETADS ) RRHO = ALPHADT / BETADS DIFFDD = 1.5e-2*0.909d0* & exp(4.6d0*exp(-0.54d0*(1.0d0/RRHO-1.0d0))) prtl_tmp = 0.15d0*RRHO elsewhere RRHO = 0.d0 DIFFDD = 0.d0 prtl_tmp = 0.d0 endwhere where (RRHO > 0.5d0) prtl_tmp = (1.85d0 - 0.85d0/RRHO)*RRHO VDC(:,:,k,1) = VDC(:,:,k,1) + DIFFDD VDC(:,:,k,2) = VDC(:,:,k,2) + prtl_tmp*DIFFDD ! if (myid==73 .and.k==4) write(*,*) 'diffdd',diffdd(20,:) end do end subroutine ddmix subroutine bldepth(trans,dbloc, dbsfc, trcr, uuu, vvv, stf, shf_qsw,& hblt, ustar, bfsfc, stable, kbl, smf, smft) ! !DESCRIPTION: ! This routine computes the ocean boundary layer depth defined as ! the shallowest depth where the bulk Richardson number is equal to ! a critical value, Ricr. ! ! NOTE: bulk richardson numbers are evaluated by computing ! differences between values at zgrid(kl) $< 0$ and surface ! reference values. currently, the reference values are equal ! to the values in the surface layer. when using higher ! vertical grid resolution, these reference values should be ! computed as the vertical averages from the surface down to ! epssfc*zgrid(kl). ! ! This routine also computes where surface forcing is stable ! or unstable (STABLE) implicit none ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1,2), intent(in) :: & trcr ! tracers at current time real*8, dimension(i0,j0,k1), intent(in) :: & uuu,vvv, &! velocities at current time dbloc, &! buoyancy difference between adjacent levels dbsfc, &! buoyancy difference between level and surface trans ! transmission (1 to 0) from surface (calculated in kppglo) real*8, dimension(i0,j0,2), intent(in) :: & stf ! surface forcing for all tracers real*8, dimension(i0,j0), intent(in) :: & shf_qsw ! short-wave forcing real*8, dimension(i0,j0,2), intent(in), optional :: & smf, &! surface momentum forcing at U points smft ! surface momentum forcing at T points ! *** either one or the other (not ! *** both) should be passed ! !OUTPUT PARAMETERS: integer, dimension(i0,j0), intent(out) :: & kbl ! index of first lvl below hbl real*8, dimension(i0,j0), intent(out) :: & hblt, &! boundary layer depth bfsfc, &! Bo+radiation absorbed to d stable, &! =1 stable forcing; =0 unstab ustar ! surface friction velocity !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & i,j,k,kkk, &! loop indices kupper, kup, kdn, ktmp, kl ! vertical level indices real*8, dimension(i0,j0) :: & vshear, &! (velocity shear re sfc)^2 sigma, &! d/hbl wm, ws, &! turb vel scale functions bo, &! surface buoyancy forcing bosol, &! radiative buoyancy forcing talpha, &! temperature expansion coeff sbeta, &! salinity expansion coeff rho1, &! density at the surface work, &! temp array zkl, &! depth at current z level b_frqncy, &! buoyancy frequency rsh_hblt, &! resolved shear contribution to HBLT (fraction) hlangm, &! Langmuir depth hekman, &! Eckman depth limit hlimit,test ! limit to mixed-layer depth ! (= min(HEKMAN,HMONOB)) real*8, dimension(i0,j0,3) :: & ri_bulk, &! Bulk Ri number at 3 lvls hmonob ! Monin-Obukhov depth limit real*8 :: & absorb_frac, &! shortwave absorption frac sqrt_arg, &! dummy sqrt argument z_upper, z_up ! upper depths for RI_BULK interpolation real*8 :: & a_co, b_co, c_co ! coefficients of the quadratic equation ! $(a_{co}z^2+b_{co}|z|+c_{co}=Ri_b) used to ! find the boundary layer depth. when ! finding the roots, c_co = c_co - Ricr real*8 :: & slope_up ! slope of the above quadratic equation ! at zup. this is used as a boundary ! condition to determine the coefficients. !----------------------------------------------------------------------- ! compute friction velocity USTAR. !----------------------------------------------------------------------- if (present(smft)) then ustar = sqrt(sqrt(smft(:,:,1)**2 + smft(:,:,2)**2)) else work = sqrt(sqrt(smf(:,:,1)**2 + smf(:,:,2)**2)) ustar = work endif !if (myid==73) write(*,*)'ustar=',ustar(20,:) !----------------------------------------------------------------------- ! compute density and expansion coefficients at surface !----------------------------------------------------------------------- work = merge(-2.0d0,trcr(:,:,1,1),trcr(:,:,1,1) < -2.0d0) call state(1,1,work,trcr(:,:,1,2), & rhofull=rho1, drhodt=talpha, drhods=sbeta) !----------------------------------------------------------------------- ! compute turbulent and radiative sfc buoyancy forcing !----------------------------------------------------------------------- !if (myid==73) write(*,*) 'stin',stf(20,4,1),stf(20,4,2),shf_qsw(20,4) !if (myid==73) write(*,*) 'stout',rho1(20,4),talpha(20,4),sbeta(20,4),grav do j=1,j0 do i=1,i0 if (rho1(i,j) /= 0.d0) then bo (i,j) = grav*(-talpha(i,j)*stf(i,j,1) - & sbeta (i,j)*stf(i,j,2))/rho1(i,j) bosol(i,j) = -grav*talpha(i,j)*shf_qsw(i,j)/rho1(i,j) else bo (i,j) = 0.d0 bosol(i,j) = 0.d0 endif end do end do ! if (myid==73) write(*,*) 'bo=',bo(20,:) ! if (myid==73) write(*,*) 'sol=',bosol(:,:) !----------------------------------------------------------------------- ! Find bulk Richardson number at every grid level until > Ricr ! max values when Ricr never satisfied are KBL = KMT and ! HBLT = -zgrid(KMT) ! ! NOTE: the reference depth is -epssfc/2.*zgrid(i,k), but the ! reference u,v,t,s values are simply the surface layer ! values and not the averaged values from 0 to 2*ref.depth, ! which is necessary for very fine grids(top layer < 2m ! thickness) ! ! ! Initialize hbl and kbl to bottomed out values ! Initialize HEKMAN and HLIMIT (= HMONOB until reset) to model bottom ! Initialize Monin Obukhov depth to value at z_up ! Set HMONOB=-zgrid(km) if unstable !----------------------------------------------------------------------- kupper = 1 kup = 2 kdn = 3 z_upper = 0.d0 z_up = zgrid(1) ri_bulk(:,:,kupper) = 0.d0 ri_bulk(:,:,kup) = 0.d0 kbl = merge( kb(:,:), 1, (kb(:,:) > 1)) hlangm = 0.d0 do kl=1,k1 zkl = -zgrid(kl) do j=1,j0 do i=1,i0 if (kl == kbl(i,j)) hblt(i,j) = zkl(i,j) end do end do enddo ! if (myid==73) write(*,*) 'testtest',hblt(20,:) if ( lcheckekmo ) then hekman = -zgrid(k1) + eps hlimit = -zgrid(k1) + eps if ( lshort_wave ) then bfsfc = bo + bosol*(1.d0-trans(:,:,1)) else bfsfc = bo endif stable = merge(1.d0, 0.d0, bfsfc >= 0.d0) bfsfc = bfsfc + stable*eps work = stable * cmonob*ustar*ustar*ustar/vonkar/bfsfc & + (stable -1.0d0)*zgrid(k1) hmonob(:,:,kup) = merge( -z_up+eps, work, work <= -z_up ) endif rsh_hblt = 0.d0 !----------------------------------------------------------------------- ! compute velocity shear squared on U-grid and use the maximum ! of the four surrounding U-grid values for the T-grid. !----------------------------------------------------------------------- do kl = 2,k1 ! loop, 1846 work = (uuu(:,:,1)-uuu(:,:,kl))**2 + & (vvv(:,:,1)-vvv(:,:,kl))**2 ZKL = -zgrid(kl) vshear(:,1) = 0.d0 do j=2,j0 vshear(1,j) = 0.d0 do i=2,i0 vshear(i,j) = max(work(i,j ), work(i-1,j ), & work(i,j-1), work(i-1,j-1)) enddo enddo !if (kl==4 .and. myid==73) write(*,*) 'trans=',trans(:,:,1) !----------------------------------------------------------------------- ! compute bfsfc= Bo + radiative contribution down to hbf * hbl ! add epsilon to BFSFC to ensure never = 0 !----------------------------------------------------------------------- if (lshort_wave) then bfsfc = bo + bosol*(1.d0-trans(:,:,kl-1)) else bfsfc = bo endif stable = merge(1.d0, 0.d0, bfsfc >= 0.d0) bfsfc = bfsfc + stable*eps !if (kl==4 .and. myid==73) write(*,*) 'bfsfc=',bfsfc(20,:) !----------------------------------------------------------------------- ! compute the Ekman and Monin Obukhov depths using above stability !----------------------------------------------------------------------- if (lcheckekmo) then fcort(1) = f(1) fcort(j0) = f(j2) do j=1,j2 fcort(j+1) = f(j) enddo do j=1,j0 do i=1,i0 if ( stable(i,j) > 0.5d0 .and. hekman(i,j) >= -zgrid(k1) ) then hekman(i,j) = max(zkl(i,j), & cekman*ustar(i,j)/(abs(fcort(j))+eps)) endif end do end do hmonob(:,:,kdn) = stable*cmonob*ustar*ustar*ustar/vonkar/bfsfc + & (stable-1.0d0)*zgrid(k1) do j=1,j0 do i=1,i0 if (hmonob(i,j,kdn) <= zkl(i,j) .and. & hmonob(i,j,kup) > -z_up) then work(i,j) = (hmonob(i,j,kdn) - hmonob(i,j,kup))/ & (z_up + zkl(i,j)) hlimit(i,j) = (hmonob(i,j,kdn) - work(i,j)*zkl(i,j))/ & (1.0d0 - work(i,j)) endif end do end do endif !----------------------------------------------------------------------- ! compute velocity scales at sigma, for hbl = -zgrid(kl) !----------------------------------------------------------------------- sigma = epssfc !if (kl==4 .and. myid==73) write(*,*) 'var',ustar,bfsfc call wscale(sigma, zkl, ustar, bfsfc, 2, wm, ws) !if (kl==4 .and. myid==73) write(*,*) 'wm=',wm(20,:) !if (kl==4 .and. myid==73) write(*,*) 'ws=',ws(20,:) !----------------------------------------------------------------------- ! compute the turbulent shear contribution to RI_BULK and store ! in WM. !----------------------------------------------------------------------- b_frqncy = sqrt( & .5d0*(dbloc(:,:,kl) + abs(dbloc(:,:,kl)) + eps2)/ & (zgrid(kl)-zgrid(kl+1)) ) wm = zkl*ws*b_frqncy* & ( (vtc/Ricr(kl))*max(2.1d0 - 200.d0*b_frqncy,concv) ) !if (kl==4 .and. myid==73) write(*,*) 'wmbq=',wm(20,:) !----------------------------------------------------------------------- ! compute bulk Richardson number at new level !----------------------------------------------------------------------- work = merge( (zgrid(1)-zgrid(kl))*dbsfc(:,:,kl), & 0.d0, kb(:,:) >= kl) ! if (kl==4 .and. myid==73) write(*,*) 'dzz=',zgrid(1)-zgrid(kl) ! if (kl==4 .and. myid==73) write(*,*) 'sfc=',dbsfc(20,:,kl) if ( linertial ) then ri_bulk(:,:,kdn) = work/(vshear+wm+ustar*bolus_sp(:,:)+eps) else ri_bulk(:,:,kdn) = work/(vshear+wm+eps) endif !if (kl==4 .and. myid==73) write(*,*) 'work=',work(20,:) !if (myid==73) write(*,*) 'bi_bulk=',kl,ri_bulk(20,:,kdn) !----------------------------------------------------------------------- ! find hbl where Rib = Ricr. if possible, use a quadratic ! interpolation. if not, linearly interpolate. the quadratic ! equation coefficients are determined using the slope and ! Ri_bulk at z_up and Ri_bulk at zgrid(kl). the slope at ! z_up is computed linearly between z_upper and z_up. ! compute Langmuir depth always !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 if ( kbl(i,j) == kb(i,j) .and. & ri_bulk(i,j,kdn) > Ricr(kl) ) then slope_up = (RI_BULK(i,j,kupper) - RI_BULK(i,j,kup))/ & (z_up - z_upper) a_co = (RI_BULK(i,j,kdn) - RI_BULK(i,j,kup) - & slope_up*(ZKL(i,j) + z_up) )/(z_up + ZKL(i,j))**2 b_co = slope_up + 2.d0 * a_co * z_up c_co = RI_BULK(i,j,kup) + & z_up*(a_co*z_up + slope_up) - Ricr(kl) sqrt_arg = b_co**2 - 4.d0*a_co*c_co if ( ( abs(b_co) > eps .and. abs(a_co)/abs(b_co) <= eps ) & .or. sqrt_arg <= 0.d0 ) then HBLT(i,j) = -z_up + (z_up + ZKL(i,j)) * & (Ricr(kl) - RI_BULK(i,j,kup))/ & (RI_BULK(i,j,kdn) - RI_BULK(i,j,kup)) else HBLT(i,j) = (-b_co + sqrt(sqrt_arg)) / (2.d0*a_co) endif kbl(i,j) = kl RSH_HBLT(i,j) = (VSHEAR(i,j)*Ricr(kl)/ & (DBSFC(i,j,kl)+eps))/HBLT(i,j) HLANGM(i,j) = USTAR(i,j) * SQRT(FSTOKES(i,j)*ZKL(i,j)/(DBSFC(i,j,kl)+eps) ) & / .9d0 endif enddo enddo !call check_val(hblt(:,:),'hblt.nc') !----------------------------------------------------------------------- ! swap klevel indices and move to next level !----------------------------------------------------------------------- ktmp = kupper kupper = kup kup = kdn kdn = ktmp z_upper = z_up z_up = zgrid(kl) end do ! end loop, 1846 !if (myid==73) write(*,*) 'hblt,bef',hblt(20,:) !----------------------------------------------------------------------- ! apply Langmuir parameterization if requested !----------------------------------------------------------------------- if ( llangmuir ) then do kl = k1,2,-1 where ( hlangm > hblt .and. & hlangm > -zgrid(kl-1) .and. & hlangm <= zkl ) hblt = hlangm kbl = kl end where enddo endif !----------------------------------------------------------------------- ! first combine Ekman and Monin-Obukhov depth limits. then apply ! these restrictions to HBLT. note that HLIMIT is set to -zgrid(km) ! in unstable forcing. !----------------------------------------------------------------------- if ( lcheckekmo ) then where ( hekman < hlimit ) hlimit = hekman do kl = 2,k1 where ( hlimit < hblt .and. & hlimit > -zgrid(kl-1) .and. & hlimit <= zkl ) hblt = hlimit kbl = kl end where enddo endif !----------------------------------------------------------------------- ! apply a Gaussian filter !----------------------------------------------------------------------- call smooth_hblt(.true., .false., hblt=hblt, kbl=kbl) !if (myid==73) write(*,*) 'hblt,aft',hblt(20,:) ! call check_val(hblt(:,:),'hbltsmo.nc') !----------------------------------------------------------------------- ! correct stability and buoyancy forcing for SW up to boundary layer !----------------------------------------------------------------------- if (lshort_wave) bfsfc = bo + bosol*(1.0d0-trans(:,:,1)) stable = merge(1.0d0, 0.d0, bfsfc >= 0.d0) bfsfc = bfsfc + stable * eps ! ensures bfsfc never=0 end subroutine bldepth subroutine wscale(sigma, hbl, ustar, bfsfc, m_or_s, wm, ws) ! !DESCRIPTION: ! Computes turbulent velocity scales. ! ! For $\zeta \geq 0, ! w_m = w_s = \kappa U^\star/(1+5\zeta)$ ! ! For $\zeta_m \leq \zeta < 0, ! w_m = \kappa U^\star (1-16\zeta)^{1\over 4}$ ! ! For $\zeta_s \leq \zeta < 0, ! w_s = \kappa U^\star (1-16\zeta)^{1\over 2}$ ! ! For $\zeta < \zeta_m, ! w_m = \kappa U^\star (a_m - c_m\zeta)^{1\over 3}$ ! ! For $\zeta < \zeta_s, ! w_s = \kappa U^\star (a_s - c_s\zeta)^{1\over 3}$ ! ! where $\kappa$ is the von Karman constant. implicit none ! !INPUT PARAMETERS: integer, intent(in) :: & m_or_s ! flag =1 for wm only, 2 for ws, 3 for both real*8, dimension(i0,j0), intent(in) :: & sigma, &! normalized depth (d/hbl) hbl, &! boundary layer depth bfsfc, &! surface buoyancy forcing ustar ! surface friction velocity ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0), intent(out) :: & wm, &! turb velocity scales: momentum ws ! turb velocity scales: tracer !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: i,j ! dummy loop indices real*8, dimension(i0,j0) :: & zeta, &! d/L or sigma*hbl/L(monin-obk) zetah ! sigma*hbl*vonkar*BFSFC or ZETA = ZETAH/USTAR**3 !----------------------------------------------------------------------- ! compute zetah and zeta - surface layer is special case !----------------------------------------------------------------------- zetah = sigma*hbl*vonkar*bfsfc zeta = zetah/(ustar**3 + eps) !----------------------------------------------------------------------- ! compute velocity scales for momentum !----------------------------------------------------------------------- if (m_or_s == 1 .or. m_or_s == 3) then do j=1,j0 do i=1,i0 if (zeta(i,j) >= 0.d0) then ! stable region wm(i,j) = vonkar*ustar(i,j)/(1.0d0 + 5.0d0*zeta(i,j)) else if (zeta(i,j) >= zeta_m) then wm(i,j) = vonkar*ustar(i,j)*(1.0d0 - 16.0d0*zeta(i,j))**0.25d0 else wm(i,j) = vonkar*(a_m*(ustar(i,j)**3)-c_m*zetah(i,j))**0.33d0 endif end do end do endif !----------------------------------------------------------------------- ! compute velocity scales for tracers !----------------------------------------------------------------------- if (m_or_s == 2 .or. m_or_s == 3) then do j=1,j0 do i=1,i0 if (zeta(i,j) >= 0.d0) then ws(i,j) = vonkar*ustar(i,j)/(1.0d0 + 5.0d0*zeta(i,j)) else if (zeta(i,j) >= zeta_s) then ws(i,j) = vonkar*ustar(i,j)*sqrt(1.0d0 - 16.0d0*zeta(i,j)) else ws(i,j) = vonkar*(a_s*(ustar(i,j)**3)-c_s*zetah(i,j))**0.33d0 endif end do end do endif end subroutine wscale subroutine smooth_hblt(overwrite_hblt, use_hmxl, & hblt, kbl, smooth_out) ! !DESCRIPTION: ! This subroutine uses a 1-1-4-1-1 Laplacian filter one time ! on HBLT or HMXL to reduce any horizontal two-grid-point noise. ! If HBLT is overwritten, KBL is adjusted after smoothing. implicit none ! !INPUT PARAMETERS: logical, intent(in) :: & overwrite_hblt, & ! if .true., HBLT is overwritten ! if .false., the result is returned in ! a dummy array use_hmxl ! if .true., smooth HMXL ! if .false., smooth HBLT ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0), optional, intent(inout) :: & hblt ! boundary layer depth integer, dimension(i0,j0), optional, intent(inout) :: & kbl ! index of first lvl below hbl ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0), optional, intent(out) :: & smooth_out ! optional output array containing the ! smoothened field if overwrite_hblt is false !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & i, j, k real*8, dimension(i0,j0) :: & work1, work2 real*8 :: & ccc, ccw, cce, ccn, ccs, & ! averaging weights, * change to ccc and ccs ztmp ! temp for level depth !----------------------------------------------------------------------- ! perform one smoothing pass since we cannot do the necessary ! boundary updates for multiple passes. !----------------------------------------------------------------------- if ( use_hmxl ) then work2 = hmxl else work2 = hblt call mpi_exch_2d_r8(work2(1,1)) endif work1 = work2 ! if(myid==73) write(*,*) 'KB=',int(kb(:,:)) ! if(myid==73) write(*,*) 'work2=',work2(20,:) do j=2,j0-1 do i=2,i0-1 if ( in(i,j,1)*kb(i,j) /= 0 ) then ccw = .125d0 cce = .125d0 ccn = .125d0 ccs = .125d0 ccc = .5d0 if ( in(i-1,j,1)*kb(i-1,j) == 0 ) then ccc = ccc + ccw ccw = 0.d0 endif if ( in(i+1,j,1)*kb(i+1,j) == 0 ) then ccc = ccc + cce cce = 0.d0 endif if ( in(i,j-1,1)*kb(i,j-1) == 0 ) then ccc = ccc + ccs ccs = 0.d0 endif if ( in(i,j+1,1)*kb(i,j+1) == 0 ) then ccc = ccc + ccn ccn = 0.d0 endif work2(i,j) = ccw * work1(i-1,j) & + cce * work1(i+1,j) & + ccs * work1(i,j-1) & + ccn * work1(i,j+1) & + ccc * work1(i,j) endif enddo enddo call mpi_exch_2d_r8(work2(:,:)) ! if (myid==73) write(*,*)'work2a',work2(20,:) ! if (myid==73) write(*,*)'zzz',-zgrid(44),-zgrid(43) do k=1,k1 do j=2,j0-1 do i=2,i0-1 ztmp = -zgrid(k) if ( k == kb(i,j) .and. work2(i,j) > ztmp ) then work2(i,j) = ztmp endif enddo enddo enddo if ( overwrite_hblt .and. .not.use_hmxl ) then hblt = work2 do k=1,k1 do j=2,j0-1 do i=2,i0-1 ztmp = -zgrid(k) if ( in(i,j,1)*kb(i,j) /= 0 .and. & ( hblt(i,j) > -zgrid(k-1) ) .and. & ( hblt(i,j) <= ztmp ) ) kbl(i,j) = k enddo enddo enddo else smooth_out = work2 endif end subroutine smooth_hblt subroutine blmix(visc, vdc, hblt, ustar, bfsfc, stable, & kbl, ghat) ! !DESCRIPTION: ! This routine computes mixing coefficients within boundary layer ! which depend on surface forcing and the magnitude and gradient ! of interior mixing below the boundary layer (matching). These ! quantities have been computed in other routines. ! ! Caution: if mixing bottoms out at hbl = -zgrid(km) then ! fictitious layer at km+1 is needed with small but finite width ! hwide(km+1). implicit none ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0,0:k0), intent(inout) :: & visc ! interior mixing coeff on input ! combined interior/bndy layer coeff output real*8, dimension(i0,j0,0:k0,2), intent(inout) :: & vdc ! diffusivity for tracer diffusion ! !INPUT PARAMETERS: integer, dimension(i0,j0), intent(in) :: & kbl ! index of first lvl below hbl real*8, dimension(i0,j0), intent(in) :: & hblt, & ! boundary layer depth bfsfc, & ! surface buoyancy forcing stable, & ! =1 stable forcing; =0 unstab ustar ! surface friction velocity ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0,k1),intent(out) :: & ghat ! non-local mixing coefficient !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer :: & k,kp1, &! dummy k level index i,j ! horizontal indices integer, dimension(i0,j0) :: & kn ! klvl closest to HBLT real*8, dimension(i0,j0,k1,3) :: & blmc ! bndy layer mixing coefs real*8, dimension(i0,j0,3) :: & gat1, &! shape function at sigma=1 dat1, &! derivative of shape function dkm1 ! bndy layer difs at kbl-1 lvl real*8, dimension(i0,j0) :: & wm,ws, &! turbulent velocity scales casea, &! =1 in case A, =0 in case B sigma, &! normalized depth (d/hbl) visch, &! viscosity at hbl difth, &! temp diffusivity at hbl difsh, &! tracer diffusivity at hbl delhat, r, dvdzup, dvdzdn, & viscp, diftp, difsp, f1, & work1,work2 !----------------------------------------------------------------------- ! compute velocity scales at hbl !----------------------------------------------------------------------- sigma = epssfc call wscale(sigma, hblt, ustar, bfsfc, 3, wm, ws) !----------------------------------------------------------------------- ! determine caseA = 0 if closer to KBL than KBL-1 ! KN is then the closest klevel to HBLT !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 k = kbl(i,j) casea(i,j) = .5d0 + SIGN(.5d0, -zgrid(k)-.5d0*hwide(k)-hblt(i,j)) enddo enddo kn = NINT(casea)*(kbl-1) + (1-NINT(casea))*kbl !----------------------------------------------------------------------- ! find the interior viscosities and derivatives at hbl by ! interpolating derivative values at vertical interfaces. compute ! matching conditions for shape function. !----------------------------------------------------------------------- f1 = stable*5.d0*bfsfc/(ustar**4+eps) do k=1,k1 do j=1,j0 do i=1,i0 if (k == kn(i,j)) then DELHAT(i,j) = .5d0*hwide(k) - zgrid(k) - HBLT(i,j) R (i,j) = 1.d0 - DELHAT(i,j) / hwide(k) DVDZUP(i,j) = (VISC(i,j,k-1) - VISC(i,j,k ))/hwide(k) DVDZDN(i,j) = (VISC(i,j,k ) - VISC(i,j,k+1))/hwide(k+1) VISCP (i,j) = .5d0*( (1.d0-R(i,j))* & (DVDZUP(i,j) + abs(DVDZUP(i,j))) + & R(i,j) * & (DVDZDN(i,j) + abs(DVDZDN(i,j))) ) DVDZUP(i,j) = (VDC(i,j,k-1,2) - VDC(i,j,k ,2))/hwide(k) DVDZDN(i,j) = (VDC(i,j,k ,2) - VDC(i,j,k+1,2))/hwide(k+1) DIFSP (i,j) = .5d0*( (1.d0-R(i,j))* & (DVDZUP(i,j) + abs(DVDZUP(i,j))) + & R(i,j) * & (DVDZDN(i,j) + abs(DVDZDN(i,j))) ) DVDZUP(i,j) = (VDC(i,j,k-1,1) - VDC(i,j,k ,1))/hwide(k) DVDZDN(i,j) = (VDC(i,j,k ,1) - VDC(i,j,k+1,1))/hwide(k+1) DIFTP (i,j) = .5d0*( (1.d0-R(i,j))* & (DVDZUP(i,j) + abs(DVDZUP(i,j))) + & R(i,j) * & (DVDZDN(i,j) + abs(DVDZDN(i,j))) ) VISCH(i,j) = VISC(i,j,k) + VISCP(i,j)*DELHAT(i,j) DIFSH(i,j) = VDC(i,j,k,2) + DIFSP(i,j)*DELHAT(i,j) DIFTH(i,j) = VDC(i,j,k,1) + DIFTP(i,j)*DELHAT(i,j) GAT1(i,j,1) = VISCH(i,j) / HBLT(i,j) /(WM(i,j)+eps) DAT1(i,j,1) = -VISCP(i,j)/(WM(i,j)+eps) + & F1(i,j)*VISCH(i,j) GAT1(i,j,2) = DIFSH(i,j) / HBLT(i,j) /(WS(i,j)+eps) DAT1(i,j,2) = -DIFSP(i,j)/(WS(i,j)+eps) + & F1(i,j)*DIFSH(i,j) GAT1(i,j,3) = DIFTH(i,j) / HBLT(i,j) /(WS(i,j)+eps) DAT1(i,j,3) = -DIFTP(i,j)/(WS(i,j)+eps) + & F1(i,j)*DIFTH(i,j) endif end do end do enddo dat1 = min(dat1,0.d0) !----------------------------------------------------------------------- ! compute the dimensionless shape functions and diffusivities ! at the grid interfaces. also compute function for non-local ! transport term (GHAT). !----------------------------------------------------------------------- do k = 1,k1 sigma = (-zgrid(k) + .5d0*hwide(k)) / hblt f1 = min(sigma,epssfc) call wscale(f1, hblt, ustar, bfsfc, 3, wm, ws) do j=1,j0 do i=1,i0 BLMC(i,j,k,1) = HBLT(i,j)*WM(i,j)*SIGMA(i,j)* & (1.d0 + SIGMA(i,j)*((SIGMA(i,j)-2.d0) + & (3.d0-2.d0*SIGMA(i,j))*GAT1(i,j,1) + & (SIGMA(i,j)-1.d0)*DAT1(i,j,1))) BLMC(i,j,k,2) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)* & (1.d0 + SIGMA(i,j)*((SIGMA(i,j)-2.d0) + & (3.d0-2.d0*SIGMA(i,j))*GAT1(i,j,2) + & (SIGMA(i,j)-1.d0)*DAT1(i,j,2))) BLMC(i,j,k,3) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)* & (1.d0 + SIGMA(i,j)*((SIGMA(i,j)-2.d0) + & (3.d0-2.d0*SIGMA(i,j))*GAT1(i,j,3) + & (SIGMA(i,j)-1.d0)*DAT1(i,j,3))) GHAT(i,j,k) = (1.d0-STABLE(i,j))* cg/(WS(i,j)*HBLT(i,j) +eps) end do end do end do !----------------------------------------------------------------------- ! find diffusivities at kbl-1 grid level !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 k = kbl(i,j) - 1 sigma(i,j) = -zgrid(k)/hblt(i,j) enddo enddo f1 = min(sigma,epssfc) call wscale(f1, hblt, ustar, bfsfc, 3, wm, ws) do j=1,j0 do i=1,i0 DKM1(i,j,1) = HBLT(i,j)*WM(i,j)*SIGMA(i,j)* & (1.d0+SIGMA(i,j)*((SIGMA(i,j)-2.d0) + & (3.d0-2.d0*SIGMA(i,j))*GAT1(i,j,1) + & (SIGMA(i,j)-1.d0)*DAT1(i,j,1))) DKM1(i,j,2) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)* & (1.d0+SIGMA(i,j)*((SIGMA(i,j)-2.d0) + & (3.d0-2.d0*SIGMA(i,j))*GAT1(i,j,2) + & (SIGMA(i,j)-1.d0)*DAT1(i,j,2))) DKM1(i,j,3) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)* & (1.d0+SIGMA(i,j)*((SIGMA(i,j)-2.d0) + & (3.d0-2.d0*SIGMA(i,j))*GAT1(i,j,3) + & (SIGMA(i,j)-1.d0)*DAT1(i,j,3))) end do end do !----------------------------------------------------------------------- ! compute the enhanced mixing !----------------------------------------------------------------------- do k=1,k1-1 where (k == (KBL - 1)) & DELHAT = (HBLT + zgrid(k))/(zgrid(k)-zgrid(k+1)) do j=1,j0 do i=1,i0 if (k == (kbl(i,j) - 1)) then BLMC(i,j,k,1) = (1.d0-DELHAT(i,j))*VISC(i,j,k) + & DELHAT(i,j) *( & (1.d0-DELHAT(i,j))**2*DKM1(i,j,1) + & DELHAT(i,j)**2*(CASEA(i,j)*VISC(i,j,k)+ & (1.d0-CASEA(i,j))*BLMC(i,j,k,1))) BLMC(i,j,k,2) = (1.d0-DELHAT(i,j))*VDC(i,j,k,2) + & DELHAT(i,j)*( & (1.d0-DELHAT(i,j))**2*DKM1(i,j,2) + & DELHAT(i,j)**2*(CASEA(i,j)*VDC(i,j,k,2)+ & (1.d0-CASEA(i,j))*BLMC(i,j,k,2))) BLMC(i,j,k,3) = (1.d0-DELHAT(i,j))*VDC(i,j,k,1) + & DELHAT(i,j) *( & (1.d0-DELHAT(i,j))**2*DKM1(i,j,3) + & DELHAT(i,j)**2*(CASEA(i,j)*VDC(i,j,k,1)+ & (1.d0-CASEA(i,j))*BLMC(i,j,k,3))) GHAT(i,j,k) = (1.d0-CASEA(i,j)) * GHAT(i,j,k) endif end do end do end do !----------------------------------------------------------------------- ! combine interior and boundary layer coefficients and nonlocal term !----------------------------------------------------------------------- do k=1,k1 do j=1,j0 do i=1,i0 if (k < kbl(i,j)) then visc(i,j,k) = blmc(i,j,k,1) vdc(i,j,k,2) = blmc(i,j,k,2) vdc(i,j,k,1) = blmc(i,j,k,3) else ghat(i,j,k) = 0.d0 endif end do end do enddo end subroutine blmix subroutine make_local_coefficients(TP,SS,PP,crho,cTP,cS,cTPTP,cSTP) ! Calculates local coefficients ! TP potential temperature [degrees, Celcius] ! SS salinity [psu, effectively ppt] ! PP pressure [bars, 10^5 Pascals] ! crho local density when TP=S=0 ! cTP gradient of local density wrt TP ! cS gradient of local density wrt S ! cTPTP .5*second derivative of local density wrt TP ! cSTP cross-derivative of local density wrt S and TP [kg/m^3/C/ppt] ! USE UNESCO_density, ONLY: density_T_S_TT_TS, temperature, & ! potential_temperature_T, potential_temperature_TT, & ! potential_temperature_TS real*8, intent(in) :: TP,SS,PP real*8, intent(out) :: cTP,cS,cTPTP,cSTP,crho real*8 dum,TT ! REAL*8 temperature,potential_temperature_T ! REAL*8 potential_temperature_TS,potential_temperature_TT ! find coeffs for rho = crho + cTP*(TP-TPr) + cS*(S-Sr) ! + cTPTP*(TP-TPr)**2 + cSTP*(S-Sr)*(TP-TPr) TT=temperature(TP,SS,PP) CALL density_T_S_TT_TS(TT,SS,PP,crho,cTP,cS,cTPTP,cSTP) dum=1.d0/potential_temperature_T(TT,SS,PP) cTP=cTP*dum cSTP = dum*(cSTP - cTP*potential_temperature_TS(PP)) cTPTP=.5d0*dum*dum*(cTPTP-cTP*potential_temperature_TT(TT,PP)) ! modify coeffs for rho = crho + cTP*TP + cS*S +cTPTP*TP**2 + cSTP*S*TP crho=crho-(cTP-cTPTP*TP)*TP-cS*SS+cSTP*SS*TP cTP=cTP-2.d0*cTPTP*TP-cSTP*SS cS=cS-cSTP*TP end subroutine make_local_coefficients !------------------------------------------------------------------------------ real*8 function temperature(TP,SS,PP) ! Calculates insitu temperature from potential temperature ! TP potential temperature [degrees Celcius] ! SS salinity [ppt] ! PP pressure [bars, 10^5 pascals] ! ! The UNESCO equation for potential_temperature is inverted iteratively. ! ! See also: potential_temperature, density REAL*8 TP,SS,PP,TT!,potential_temperature INTEGER*2 iterate TT=TP temperature=TT ! first `guess' ! under relax to obtain temperature from potential temperatuer DO iterate=1,100 TT=TT+.999d0*(TP-potential_temperature(TT,SS,PP)) ! T=T+(TP-potential_temperature(T,S,P)) IF (DABS(temperature-TT)<1.D-5) THEN ! PRINT *,'number of iterations = ',iterate EXIT ENDIF temperature=TT ENDDO IF(iterate>=100) THEN PRINT *,'ERROR: temperature did not converge in 100 iterations' PRINT *,'Program stopped by FUNCTION temperature in indata' CALL MPI_FINALIZE(IERR) STOP 0 ENDIF end function temperature !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ real*8 function potential_temperature(TT,SS,PP) ! Calculates the potential temperature from ! T temperature [degrees Celcius] ! S salinity in [parts per thousand (ppt)] ! P pressure in [bars (10^5 Pascals)] ! EXAMPLE: ! t=potential_temperature(T=10,S=25,P=1000) ! gives t=8.46785160000000 ! References: ! Gill, A. E. (1981) Atmosphere-Ocean Dynamics. Academic Press ! Bryden, H. L. (1973) New polynomials for thermal expansion, ! adiabatic temperature gradient and geopotential ! temperature gradient of sea water. Deep-Sea Res., 20, 401-408. REAL*8 TT,SS,PP potential_temperature = & TT-PP*(3.6504D-4 + TT*(8.3198D-5 - TT*(5.4065D-7-TT*4.0274D-9))) & -PP*(SS-35.d0)*(1.7439D-5 - TT*2.9778D-7) & -PP*PP*(8.9309D-7 - TT*(3.1628D-8 - TT*2.1987D-10)) & +4.1057D-9*(SS-35.d0)*PP*PP & -PP*PP*PP*(-1.6056D-10 + TT*5.0484D-12) end function potential_temperature !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ real*8 function potential_temperature_T(TT,SS,PP) ! Calculates the derivative of potential temperature wrt insitu temperature T ! TT temperature in degrees Celcius ! SS salinity in parts per thousand ! PP pressure in bars (10^5 Pascals) ! ! Note: the derivative of T wrt potential temperature is ! 1/potential_temperature_T ! ! References: ! Gill, A. E. (1981) Atmosphere-Ocean Dynamics. Academic Press ! Bryden, H. L. (1973) New polynomials for thermal expansion, ! adiabatic temperature gradient and geopotential ! temperature gradient of sea water. Deep-Sea Res., 20, 401-408. REAL*8 TT,SS,PP potential_temperature_T= & 1.+PP*(-8.3198e-5+TT*(1.0813e-6-TT*1.20822e-8) & +(SS-35.)*2.9778e-7+PP*((3.1628e-8-TT*4.3974e-10)-PP*5.0484e-12)) end function potential_temperature_T !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ real*8 function potential_temperature_TT(TT,PP) ! Calculates the 2nd derivative of potential temperature wrt temperature T ! TT temperature in degrees Celcius ! PP pressure in bars (10^5 Pascals) ! ! References: ! Gill, A. E. (1981) Atmosphere-Ocean Dynamics. Academic Press ! Bryden, H. L. (1973) New polynomials for thermal expansion, ! adiabatic temperature gradient and geopotential ! temperature gradient of sea water. Deep-Sea Res., 20, 401-408. REAL*8 TT,PP potential_temperature_TT = & PP*( 1.0813e-6 - TT*2.41644e-8 - PP*4.3974e-10 ) end function potential_temperature_TT !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ real*8 function potential_temperature_TS(PP) ! Calculates the second derivative of potential temperature wrt ! insitu temperature T and salinity S ! P pressure in bars (10^5 Pascals) ! ! References: ! Gill, A. E. (1981) Atmosphere-Ocean Dynamics. Academic Press ! Bryden, H. L. (1973) New polynomials for thermal expansion, ! adiabatic temperature gradient and geopotential ! temperature gradient of sea water. Deep-Sea Res., 20, 401-408. REAL*8 PP potential_temperature_TS = PP*2.9778e-7 end function potential_temperature_TS !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ subroutine density_T_S_TT_TS(TT,SS,PP,rho,drho_dT,drho_dS, & drho_dT_dT,drho_dT_dS) ! Uses equation (A3.1) of Gill 1982 Atmosphere Ocean Dynamics to ! calculate the density of fresh water and then uses (A3.2) to ! calculate the density of saline water at one standard atmosphere. ! Pressure effects then accounted for using (A3.3). ! ! The derivatives of density with respect to temperature T and ! salinity S are also calculated. The second derivative wrt T ! and the TS cross-derivative are calculated. ! ! INPUT VARIABLES: ! TT temperature, [degrees Celcius] ! SS salinity, [psu, which effectively equals ppt] ! PP pressure, [bars, 10^5 Pascals] ! ! OUTPUT VARIABLES: ! rho density at T,S,P [kg/m^3] ! drho_dT derivative of density wrt T [kg/m^3/C] ! drho_dS derivative of density wrt S [kg/m^3/ppt] ! ! OPTIONAL OUTPUT VARIABLES: ! drho_dT_dT 2nd derivative of density wrt T [kg/m^3/C/C] ! drho_dT_dS cross derivative of density wrt T and S [kg/m^3/C/ppt] ! ! test values are rho(T=5, S=0, P=0) = 999.96675 - 1000 ! rho(T=5, S=35,P=0) = 1027.67547 - 1000 ! rho(T=25, S=35, P=1000) = 1062.53817 - 1000 ! ! This algorithm is modified so it works well at single precision. ! This is done by extracting an offset density of 1000 kg/m^3 ! USE global_parameters, ONLY: ks REAL*8 TT,SS,PP,rho,drho_dT,drho_dS,drho_dT_dT,drho_dT_dS REAL*8 K,dum, dKdT, dKdS, r,dKdTdT, rdT, rdS,dKdSdT ! calculate the density of fresh water at temperature TT rho = -0.1574060 + (6.793952E-2 - (9.095290E-3-(1.001685E-4 & - (1.120083E-6 - 6.536332E-9*TT)*TT)*TT)*TT)*TT ! calculate the derivative wrt T of density of fresh water at temperature T drho_dT = 6.793952E-2 - (1.819058E-2 - (3.005055E-4 & -(4.480332E-6 -3.268166E-8*TT)*TT)*TT)*TT ! calculate the 2nd derivative wrt T of density of fresh water at temperature T !CW IF(PRESENT(drho_dT_dT)) drho_dT_dT = - (1.819058E-2 - (6.01011E-4 & - (1.3440996E-5 - 1.3072664E-7*TT)*TT)*TT) dum=DBLE(SQRT(SS)) ! calculate the density of saline water at 1 atmosphere rho = rho + SS*(0.824493 +4.8314E-4*SS - (4.0899E-3 - (7.6438E-5 & - (8.2467E-7 - 5.3875E-9*TT)*TT)*TT)*TT) & + SS*dum*(-5.72466E-3 + (1.0227E-4 - 1.6546E-6*TT)*TT) ! calculate the derivative wrt T of density of saline water at 1 atmosphere drho_dT = drho_dT & +SS*(-4.0899E-3 +(1.52876E-4-(2.47401E-6 - 2.155E-8*TT)*TT)*TT) & +SS*dum*(1.0227E-4-3.3092E-6*TT) ! calculate the 2nd derivative wrt T of density of saline water at 1 atmosphere !CW IF(PRESENT(drho_dT_dT)) drho_dT_dT = drho_dT_dT & +SS*(1.52876E-4 - (4.94802E-6-6.465E-8*TT)*TT)-SS*dum*3.3092E-6 ! calculate the derivative wrt S of density of saline water at 1 atmosphere drho_dS = 0.824493-(4.0899E-3 - (7.6438E-5 & -(8.2467E-7-5.3875E-9*TT)*TT)*TT)*TT & +dum*(-8.58699E-3+(1.53405E-4-2.4819E-6*TT)*TT)+9.6628E-4*SS ! calculate the density cross-derivative wrt S and T of saline water at 1 atmos !CW IF(PRESENT(drho_dT_dS)) drho_dT_dS = - (4.0899E-3 - (1.52876E-4 & -(2.47401E-6-2.155E-8*TT)*TT)*TT)+dum*(1.53405E-4-4.9638E-6*TT) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%% PRESSURE EFFECTS via secant bulk modulus %%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! calculate the secant bulk modulus of pure water K = 19652.21 + TT*(148.4206 - TT*(2.327105 - TT*(1.360477E-2 & - TT*5.155288E-5))) ! calculate the derivative wrt T of the secant bulk modulus of pure water dKdT = 148.4206 - TT*(4.65421 - TT*(4.081431E-2 -TT*2.0621152E-4)) ! calculate the 2nd derivative wrt T of the secant bulk modulus of pure water dKdTdT = - (4.65421 - TT*(8.162862E-2 - TT*6.1863456E-4)) ! calculate the bulk modulus at one atmosphere K = K+SS*(54.6746 - TT*(.603459 -TT*(1.09987E-2 - TT*6.1670E-5))) & +SS*dum*(7.944E-2 + TT*(1.6483E-2 - TT*5.3009E-4)) ! calculate derivative wrt T of the bulk modulus at one atmosphere dKdT = dKdT +SS*(- .603459 + TT*(2.19974E-2 - TT*1.8501E-4)) & +SS*dum*(1.6483E-2 - 1.06018E-3*TT) ! calculate 2nd derivative wrt T of the bulk modulus at one atmosphere !CW IF(PRESENT(drho_dT_dT)) dKdTdT = dKdTdT +SS*(2.19974E-2 - TT*3.7002E-4 - dum*1.06018E-3) ! calculate the derivative wrt S of the bulk modulus at one atmosphere dKdS = 54.6746 - TT*(.603459 - TT*(1.09987E-2 - TT*6.1670E-5)) & +dum*(1.1916E-1 + TT*(2.47245E-2 - TT*7.95135E-4)) ! calculate the cross derivative wrt S and T ! of the bulk modulus at one atmosphere !CW IF(PRESENT(drho_dT_dS)) dKdSdT = - (.603459 - TT*(2.19974E-2 - TT*1.8501E-4)) & + dum*(2.47245E-2 - TT*1.59027E-3) ! correct the bulk modulus for pressure p K =K+PP*((3.239908+TT*(1.43713E-3+TT*(1.16092E-4-TT*5.77905E-7))) & + SS*(2.2838E-3 - TT*(1.0981E-5+TT*1.6078E-6))+1.91075E-4*SS*dum & + PP*( (8.50935E-5 - TT*(6.12293E-6 - TT*5.2787E-8)) & + SS*(-9.9348E-7 + TT*(2.0816E-8 + TT*9.1697E-10)))) ! correct the derivative wrt T of the bulk modulus for pressure p dKdT=dKdT + PP*((1.43713E-3 + TT*(2.32184E-4 - TT*1.733715E-6)) & + SS*(- 1.0981E-5 - 3.2156E-6*TT) & + PP*( (- 6.12293E-6 + 1.05574E-7*TT) & + SS*(2.0816E-8 + 1.83394E-9*TT))) ! correct the 2nd derivative wrt T of the bulk modulus for pressure p !CW IF (PRESENT(drho_dT_dT)) dKdTdT = dKdTdT + PP*( (2.32184E-4 - TT*3.46743E-6) & - SS*3.2156E-6 + PP*(1.05574E-7+ SS*1.83394E-9)) ! correct the derivative wrt S of the bulk modulus for pressure p dKdS = dKdS + PP*((2.2838E-3 - TT*(1.0981E-5 + TT*1.6078E-6)) & +2.866125E-4*dum+PP*(-9.9348E-7+TT*(2.0816E-8 + TT*9.1697E-10))) ! correct the cross derivative wrt S and T of the bulk modulus for pressure p !CW IF(PRESENT(drho_dT_dS)) dKdSdT = dKdSdT - PP*(1.0981E-5 + TT*3.2156E-6 & - PP*(2.0816E-8 + TT*1.83394E-9)) ! correct the density for pressure effects r=rho rdT=drho_dT rdS=drho_dS dum=1.d0/(K-PP) rho = (K*rho+1000*PP)*dum drho_dT = (K*drho_dT + (r - rho)*dKdT)*dum drho_dS = (K*drho_dS + (r - rho)*dKdS)*dum !CW IF(PRESENT(drho_dT_dT)) drho_dT_dT = dum*(r*dKdTdT+2*rdT*dKdT+K*drho_dT_dT & + rho*(2*dum*dKdT*dKdT-dKdTdT)- 2*dum*(r*dKdT+K*rdT)*dKdT) !CW IF(PRESENT(drho_dT_dS)) drho_dT_dS = dum*( dKdS*rdT+K*drho_dT_dS+rdS*dKdT & +r*dKdSdT-rho*dKdSdT-drho_dS*dKdT -drho_dT*dKdS) end subroutine density_T_S_TT_TS !----------------------------------------------------------------------------------- !----------------------------------------------------------------------------------- real*8 function local_density(TP,SS,crho,cTP,cS,cTPTP,cSTP) ! Calculates the density from local coefficients. ! TP potential temperature [degrees, Celcius] ! SS salinity [psu, effectively ppt] ! crho local density when TP=S=0 ! cTP gradient of local density wrt TP [kg/m^3/C] ! cS gradient of local density wrt S [kg/m^3/ppt] ! cTPTP .5*second derivative of local density wrt TP [kg/m^3/C/C] ! cSTP cross-derivative of local density wrt S and TP [kg/m^3/C/ppt] ! Operation count = 8 real*8 tp,ss,crho,cTP,cS,cTPTP,cSTP local_density=crho+(cTP+cTPTP*tp+cSTP*ss)*tp+cS*ss end function local_density !----------------------------------------------------------------------------------- !----------------------------------------------------------------------------------- subroutine ice_formation use time_management, only: ice_ts, nsteps_run use ice, only: AQICE, QFLUX, QICE use blocks, only: nx_block, ny_block !This subroutine compute ice formation and the heat flux associated with !ice formation. This heat flux is sent to the ice model via the flux coupler. !Merged and modified from pop2 ice module. implicit none ! real*8, dimension(i0,j0), save :: & ! qice, &!tot column cooling from ice form (in C*cm) ! aqice, &!sum of accumulated ice heat flux since tlast ! qflux !internal ocn heat flux due to ice formation real*8, dimension(i0,j0) :: & fw_freeze, & !water flux at T points due to frazil ice formation tfrz, & !freezing temperature of water in deg C potice, & !potential amount of ice formation work1 !working variable integer, save :: & ice_flag, &!time flag id for ice formation kmxice !lowest level from which to integrate ice formation real*8, save :: & ref_val,&!tracer reference value salice, &!sea ice salinity in msu salref, &!ocean ref salinity in msu cp_over_lhfusion !cp_sw/latent_heat_fusion = !rho_sw*cp_sw / (latent_heat_fusion*rho_fw) integer :: i,j,k !WRITE(*,*)'[physics]ice_ts=',ice_ts,nsteps_run !======================== if ( ice_ts ) then !======================== !=================================================================================== ! Initialization (first call) !=================================================================================== if (liceform_ini) then cp_over_lhfusion = (4.1d0/3.996d0) * 3.996e7 & / (3.34e9 * 1.0d0) kmxice = 1 ! ice_freq_opt = 'never' ! ice_freq = 100000 salice = sea_ice_salinity*ppt_to_salt salref = ocn_ref_salinity*ppt_to_salt ! tlast_ice = 0.0 ! qice = 0.0 ! aqice = 0.0 ! qflux = 0.0 ! fw_freeze = 0.0 liceform_ini=.false. endif QICE(:,:,1) = 0.0d0 potice = 0.0d0 p0lf(:,:) = p0(:,:) !=================================================================================== ! ice formation !=================================================================================== ! compute frazil ice formation for sub-surface layers. if ice ! forms in lower layers but layers above are warm - the heat is ! used to melt the ice. the ice formation occurs at salinity, Si. ! this volume is replaced with an equal volume at the salinity of ! the layer above. the total ice heat flux is accumulated. ! ! WARNING: unless a monotone advection scheme is in place, ! advective errors could lead to temps that are far below freezing ! in some locations and this scheme will form lots of ice. ! ice formation should be limited to the top layer (kmxice=1) ! if the advection scheme is not monotone. ! ! modified from pop2 ice.F90, partial_bottom_cells == .False. ! lfw_as_salt_flx == .True. ! sfc_layer_type == sfc_layer_varthick !----------------------------------------------------------------------- !------------------------------------------------------------------------ ! assign working variable and change the salinity unit to (g/g) !------------------------------------------------------------------------ trcr(:,:,:,1) = t2 trcr(:,:,:,2) = s2*ppt_to_salt do k=kmxice,2,-1 !*** !*** potice is the potential amount of ice formation !*** (potice>0) or melting (potice<0) in layer k !*** call tfreez(tfrz,trcr(:,:,k,2)) where ( k <= int(kb(:,:)) ) & potice = (tfrz - trcr(:,:,k,1))/odz(k) !*** !*** if potice < 0, use the heat to melt any ice !*** from lower layers !*** if potice > 0, keep on freezing (QICE < 0) !*** ! potice = max(potice,qice(:,:)) potice = max(potice,QICE(2:nx_block-1,2:ny_block-1,1)) !*** !*** adjust tracer values based on freeze/melt (tfrz or ori-val) !*** trcr(:,:,k,1) = trcr(:,:,k,1) + potice*odz(k) ref_val = salref - salice if (ref_val /= 0.0) then trcr(:,:,k,2) = trcr(:,:,k,2)+ref_val*potice*cp_over_lhfusion*odz(k) end if !*** accumulate freezing potential ! qice(:,:) = qice(:,:) - potice QICE(2:nx_block-1,2:ny_block-1,1) = QICE(2:nx_block-1,2:ny_block-1,1) & - potice end do ! k loop !----------------------------------------------------------------------- ! now repeat the above algorithm for the surface layer. when fresh ! water flux formulation is used, the surface layer does not get ! any salt from other layers. instead, its volume changes. !----------------------------------------------------------------------- k = 1 call tfreez(tfrz,trcr(:,:,k,2)) work1 = 1.d0/odz(k) ! work1 = work1 + PSURF(:,:,newtime,bid)/grav + eps2 work1 = work1 + p0lf(:,:)/grav + eps2 where ( k <= int(kb(:,:)) ) potice = (tfrz - trcr(:,:,k,1))*work1 end where ! potice = max(potice, qice(:,:)) potice = max(potice,QICE(2:nx_block-1,2:ny_block-1,1)) trcr(:,:,k,1) = trcr(:,:,k,1) + potice/work1 ref_val = salref - salice if (ref_val /= 0.0) then trcr(:,:,k,2) = trcr(:,:,k,2) + ref_val*potice*cp_over_lhfusion/work1 end if ! qice(:,:) = qice(:,:) - potice QICE(2:nx_block-1,2:ny_block-1,1) = QICE(2:nx_block-1,2:ny_block-1,1) & - potice !----------------------------------------------------------------------- ! let any residual heat in the upper layer melt previously formed ice !----------------------------------------------------------------------- ! aqice(:,:) = aqice(:,:) + time_weight*qice(:,:) AQICE(2:nx_block-1,2:ny_block-1,1) = AQICE(2:nx_block-1,2:ny_block-1,1) & + time_weight*QICE(2:nx_block-1,2:ny_block-1,1) !----------------------------------------------------------------------- ! recalculate freezing potential based on adjusted T. ! only interested in melt potential now (POTICE < 0) - use this ! melt to offset any accumulated freezing (AQICE < 0) and ! adjust T and S to reflect this melting. when freshwater flux ! formulation, compute the associated freshwater flux instead of ! adjusting S. !----------------------------------------------------------------------- call tfreez(tfrz,trcr(:,:,k,2)) where (k <= int(kb(:,:)) ) potice = (tfrz - trcr(:,:,k,1)) * work1 endwhere ! potice = max(potice, aqice(:,:)) potice = max(potice,AQICE(2:nx_block-1,2:ny_block-1,1)) trcr(:,:,k,1) = trcr(:,:,k,1) + potice/work1 ref_val = salref - salice if (ref_val /= 0.0) then trcr(:,:,k,2) = trcr(:,:,k,2) + ref_val*potice*cp_over_lhfusion/work1 end if ! aqice(:,:) = aqice(:,:) - time_weight*potice AQICE(2:nx_block-1,2:ny_block-1,1) = AQICE(2:nx_block-1,2:ny_block-1,1) & - time_weight*potice !--------------------------------------------------------------------------- ! update T2, S2 and change the salinity unit to (g/kg) !--------------------------------------------------------------------------- t2 = trcr(:,:,:,1) s2 = trcr(:,:,:,2)/ppt_to_salt endif ! time to do ice !=================================================================================== ! sent fluxes to coupler ( if needed !=================================================================================== ! if (mod(itff,itfday)==0) then ! if (lcoupling ) then call ice_flx_to_coupler end subroutine ice_formation subroutine ice_flx_to_coupler use ice, only: AQICE, QFLUX, QICE, tlast_ice use blocks, only: nx_block, ny_block use time_management, only: nsteps_run implicit none !----------------------------------------------------------------------- ! ! This subroutine sets up the ice formation / melting potential ! heat fluxes to be sent to the coupler. ice formation heat flux ! is accumulated for time averaging. ! !----------------------------------------------------------------------- integer :: k real*8, dimension(i0,j0) :: & tfrz, &! freezing temp of water work1, work2 ! work arrays integer :: i,j !----------------------------------------------------------------------- ! ! compute the first layer thickness ! !----------------------------------------------------------------------- trcr(:,:,:,1) = t2 trcr(:,:,:,2) = s2*ppt_to_salt k = 1 call tfreez(tfrz,trcr(:,:,k,2)) work1 = 1.d0/odz(k) ! WORK1 = WORK1 + PSURF(:,:,curtime,bid)/grav + eps2 work1 = work1 + p0lf(:,:)/grav + eps2 !----------------------------------------------------------------------- ! ! first compute the melt potential ! !----------------------------------------------------------------------- work2 = 0.d0 where ( k <= int(kb(:,:)) ) work2 = (tfrz - trcr(:,:,k,1)) * work1 endwhere !----------------------------------------------------------------------- ! ! adjust ice formation amount when mixing step is tavg ! !----------------------------------------------------------------------- AQICE(:,:,1) = 0.5d0 * AQICE(:,:,1) !----------------------------------------------------------------------- ! ! merge the ice formation and melt potential fluxes and compute ! !----------------------------------------------------------------------- where ( AQICE(2:nx_block-1,2:ny_block-1,1) < 0.d0 ) QICE(2:nx_block-1,2:ny_block-1,1) = -AQICE(2:nx_block-1,2:ny_block-1,1) elsewhere QICE(2:nx_block-1,2:ny_block-1,1) = work2 endwhere if (tlast_ice == 0.d0) then QFLUX(:,:,1) = 0.d0 else QFLUX(:,:,1) = QICE(:,:,1)/tlast_ice/hflux_factor endif end subroutine ice_flx_to_coupler subroutine tfreez(TFRZ,SALT) ! !DESCRIPTION: ! This function computes the freezing point of salt water. real*8, dimension(i0,j0), intent(in) :: & SALT !salinity in model units (g/g) real*8, dimension(i0,j0), intent(out) :: & TFRZ !freezing temperature of water in deg C !TFRZ = -0.0544d0*SALT*salt_to_ppt TFRZ = -1.8d0 end subroutine tfreez subroutine hv_hmix_tracer_type_gm(k,vdc) implicit none integer, intent(in) :: k real*8, dimension(i0,j0,0:k0,2), intent(inout) :: vdc integer :: i,j,l real*8 :: emax,dtmp,hm if ( hmix_tracer_type_gm > 0.5) then l=k+1 !TS emax=min(.5d0/(dt*odzw(l)**2),3000.d0) emax=.5d0/(dt*odzw(l)**2) !TSb1m0 if(k.le.3)then !TSb1m0 hm=.01d0 !TSb1m0 else !TSb1m0 hm=1.d0 !TSb1m0 endif !TS emax=MIN(1200.d0,emax) !emax=500.d0 do j=1,j2 do i=1,i2 ! if(k .ge. kbl(i+1,j+1))then ! emax=.5d0/(dt*odzw(l)**2) ! emax=MIN(6000.d0,emax) ! else ! emax=600.d0 ! endif dtmp=odzw(l)*iw(i+1,j+1,l) !! limit the value not to over 900 or .5/(DT*ODZW(L)**2) hv(i,j,k)=dtmp*MAX(MIN(emax,vdc(i+1,j+1,k,1)+add(i,j,k)),hbk(k)*5.d0) vdc(i+1,j+1,k,1)=MAX(MIN(emax,vdc(i+1,j+1,k,1)+add(i,j,k)),hbk(k)*5.d0) vdc(i+1,j+1,k,2)=MAX(MIN(emax,vdc(i+1,j+1,k,2)+add(i,j,k)),hbk(k)*.05d0) !TSb1m0 vdc(i+1,j+1,k,2)=MAX(MIN(emax,vdc(i+1,j+1,k,2)*hm+add(i,j,k)),hbk(k)*.05d0) vdct(i,j,k)=vdc(i+1,j+1,k,1) vdcs(i,j,k)=vdc(i+1,j+1,k,2) end do end do endif end subroutine hv_hmix_tracer_type_gm subroutine set_surface_forcing use time_management use forcing_coupled, only: tday00_interval_beg, interval_cum_dayfrac,& QSW_COSZ_WGHT_NORM, QSW_COSZ_WGHT, compute_cosz integer :: index_qsw real*8 :: cosz_day index_qsw = mod(nsteps_this_interval,nsteps_per_interval) + 1 cosz_day = tday00_interval_beg + interval_cum_dayfrac(index_qsw-1) & - interval_cum_dayfrac(nsteps_per_interval) call compute_cosz(cosz_day, 1, QSW_COSZ_WGHT(:,:,1)) where (QSW_COSZ_WGHT_NORM(:,:,1) > 0.0d0) QSW_COSZ_WGHT(:,:,1) = QSW_COSZ_WGHT(:,:,1) & * QSW_COSZ_WGHT_NORM(:,:,1) elsewhere QSW_COSZ_WGHT(:,:,1) = 1.0d0 endwhere ! SHF_QSW(:,:) = QSW_COSZ_WGHT(2:nx_block-1,2:ny_block-1,1)*SHF_QSW(:,:) end subroutine set_surface_forcing subroutine precip_adjustment_init use time_management, only: iyear,imonth,iday ! real*8, intent(inout) :: saltotal_tendency integer :: ymd,ITFF real*8 :: saltotal_initial ymd = iyear*10000 + imonth*100 + iday ITFF = ITF - IT0 !!!!! Part of writing saltotal_initial if((lrstrt==1.and.ITF.eq.1).or.(lrstrt==0.and.ITF.eq.181).or. & (mod(ymd-101,10000)==0 .and. mod(ITFF,ITFDAY)==1))THEN !! for a start-up case of timcom (ITF==1) OR !! the 1st time step of every Jan 1st TOTALS=SUM(S1(2:I1,2:J1,:)*VOLUME(2:I1,2:J1,:)) DTMP=TOTALS CALL MPI_ALLREDUCE(DTMP,TOTALS,1,MPI_REAL8,MPI_SUM,M_CART,IERR) saltotal_initial=TOTALS ! global total salt is saved in saltotal_initial if (myid==0) then ! use one cpu to write write(*,*) 'write sal initial at the beginning',ymd,'ITF=',ITF open(9093,file='global_ssh_saltotal_initial',form='unformatted') write(9093) saltotal_initial close(9093) WRITE(*,*)'saltotal_init0:YMD=',ymd,'ITF=',ITF,'ITFF=',ITFF WRITE(*,*)'saltotal_init0 (init ):',saltotal_initial,saltotal_initial/TOTALVOL WRITE(*,*)'saltotal_init0 (tendency):',saltotal_tendency endif endif !! ITF.eq.1 end !!!!! Part end of writing saltotal_initial end subroutine precip_adjustment_init subroutine precip_adjustment use time_management, only: iyear,imonth,iday ! real*8, intent(inout) :: saltotal_tendency integer :: ymd,ITFF REAL*8 :: saltotal_initial,saltotal_final ymd = iyear*10000 + imonth*100 + iday ITFF = ITF - IT0 !!!!!!!!!!!!!!! Section of salt adjustment if (ladjust_salinity_budget) then !!!!! Part of reading saltotal_tendency (for each submission) ! if(ITFF.eq.1 .and. ymd>11231) then if(mod(ymd-101,10000)==0 .and. ymd>11231 .and. mod(ITFF,ITFDAY)==1) then !! READ saltotal_tendency for the first time-step of a resubmission !after 1 year, !! else saltotal_tendency = 0 open(9093,file='global_saltotal_tendency',form='unformatted') read(9093) saltotal_tendency close(9093) if(myid.eq.0)then WRITE(*,*)'saltotal_init1:YMD=',ymd,'ITF=',ITF,'ITFF=',ITFF WRITE(*,*)'saltotal_init1 (init ):',saltotal_initial,saltotal_initial/TOTALVOL WRITE(*,*)'saltotal_init1 (tendency):',saltotal_tendency endif endif !! if(ladjust_salinity_budget !!!!! Part end of reading saltotal_tendency (for each submission) !!!!! Part of (1) read saltotal_initial, (2) cal. & (3) write !saltotal_tendency, on Dec. 31st if (mod(ymd-1231,10000)==0 .and. mod(ITFF,ITFDAY)==0) then !! Last day of every year !ck if (Nd==INT(DMONTH(Nm)) .and. mod(ITFF,ITFDAY)==0 !! Last day & timestep of every month !ck .and. ladjust_salinity_budget) then !! mod(ITFF,ITFDAY)==0 is for the last time-step of a day if(myid==0) write(*,*) 'read sal initial in ', ymd,'ITF=',ITF open(9093,file='global_ssh_saltotal_initial',form='unformatted') read(9093) saltotal_initial !! (1) close(9093) if(myid==0)then WRITE(*,*)'saltotal_final0:YMD=',ymd,'ITF=',ITF,'ITFF=',ITFF WRITE(*,*)'saltotal_final0 (init ):',saltotal_initial,saltotal_initial/TOTALVOL WRITE(*,*)'saltotal_final0 (tendency):',saltotal_tendency endif IF(ymd.ge.21231)then TOTALS=SUM(S1(2:I1,2:J1,:)*VOLUME(2:I1,2:J1,:)) DTMP=TOTALS CALL MPI_ALLREDUCE(DTMP,TOTALS,1,MPI_REAL8,MPI_SUM,M_CART,IERR) saltotal_final=TOTALS saltotal_tendency=saltotal_tendency+& (saltotal_final-saltotal_initial)/TOTALVOL/dble(365*DAODT) !! (2) endif ! 1 (saltotal_final-saltotal_initial)/TOTALVOL/dble(DMONTH(Nm)*DAODT) if(myid==0)THEN WRITE(*,*)'saltotal_final1 (init ):',saltotal_final,saltotal_final/TOTALVOL WRITE(*,*)'saltotal_final1 (tendency):',saltotal_tendency open(9093,file='global_saltotal_tendency',form='unformatted') write(9093) saltotal_tendency !! (3) close(9093) endif endif !! Last day of every year !!!!! Part end of (1)(2)(3) !CTS impose the correction DO K=1,K1 DO J=2,J1 DO I=2,I1 DTMP=IN(I,J,K)*saltotal_tendency!/(ITFDAY*365) S1(I,J,K)=S1(I,J,K)-DTMP S2(I,J,K)=S2(I,J,K)-DTMP SLF(I,J,K)=SLF(I,J,K)-DTMP enddo enddo enddo endif !! ladjust_salinity_budget end !!!!!!!!!!!!!!! Section of salt adjustment end end subroutine precip_adjustment end module physics