c------------------------------------------------------------
      subroutine tau_special(nstep,nstart,delt,npt,taux,tauy)
c------------------------------------------------------------
      dimension taux(npt),tauy(npt)

      PARAMETER(JPT=31)
      common/mlevy/ztaux,ztauy,ztot,zsol
      REAL ZTAUX(JPT),ZTAUY(JPT),ZTOT(JPT),ZSOL(JPT)

      this_time = 1. + (nstep-nstart)*delt
      i_time = int(this_time)
      al = this_time - i_time
      this_taux = ((1.-al)*ztaux(i_time) + al*ztaux(i_time+1))/1030.
      this_tauy = ((1.-al)*ztauy(i_time) + al*ztauy(i_time+1))/1030.
      do i = 1, npt 
         taux(i) = this_taux
         tauy(i) = this_tauy
      enddo

      return
      end

c------------------------------------------------------------
      subroutine qforc_special(nstep,nstart,delt,npt,q,qr,qsr)
c------------------------------------------------------------
c  q+ qr = total heat budget
c  qr is the penetrative part of the solar radiation (i. e. qr=0 is solar penetration
c       is not allowed)
c q and qr are in dT/dt units (and not in W/m2)
c qsr is only used for the biology : it is the total solar radiation, in W/m2
c if use_ice, qsr is computed within amlice_flux
c
      dimension q(npt),qr(npt),qsr(npt)

      PARAMETER(JPT=31)
      common/mlevy/ztaux,ztauy,ztot,zsol
      REAL ZTAUX(JPT),ZTAUY(JPT),ZTOT(JPT),ZSOL(JPT)

      common /new_hfxevp/ QCON, rlx_time, solr_gamma, TATM, SATM
     *                    ,trans_coef_sst,trans_coef_sss,coef_precip

      qcon_gam = solr_gamma/QCON
      qcon_inv = 1/QCON
      this_time = 1. + (nstep-nstart)*delt
      i_time = int(this_time)
      al = this_time - i_time
      this_qtot = (1.-al)*ztot(i_time) + al*ztot(i_time+1)
      this_qsol = (1.-al)*zsol(i_time) + al*zsol(i_time+1)
      do i = 1, npt 
         qr(i)= this_qsol * qcon_gam
         q(i) = this_qtot * qcon_inv - this_qsol * qcon_gam
      enddo
      do i = 1, npt 
         qsr(i) = qr(i)/ qcon_gam
      enddo
      
      return
      end

c------------------------------------------------------------
      subroutine misc_out (tenso)
c------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_nao.h'

      if (iep_key.gt.0) then
         write(IEP_OUT,303)coef_bud,budget_eb,budget_pr,
     *        budget_sp,budget_ri,budget_re,bud
         call flush (IEP_OUT)
      endif
      if (inao_key.gt.0) then
         write(INAO_OUT,*)tenso,tscl_anom,tscl_a
         call flush (INAO_OUT)
      endif

  303 format(2x,f8.6,4x,6e12.3)

      return
      end

c------------------------------------------------------------
      subroutine bud_reset
c------------------------------------------------------------
      include 'comm_new.h'

      budget_eb = 0.
      budget_pr = 0.
      budget_sp = 0.
      budget_ri = 0.
      budget_re = 0.

      return
      end

c------------------------------------------------------------
      subroutine set_ramp(delt,nstep,w_day)
c------------------------------------------------------------
      common /test_time/ t_day,t_ramp

      w_day = delt*nstep
      t_day  = w_day - int(w_day)
      t_ramp = min(1.,w_day/20.)

      return
      end

c------------------------------------------------------------
      subroutine force_bio(npt,nz,nstep,flor,tmonth)
c------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_bio.h'
      include 'comm_data.h'
      include 'comm_tracer.h'

      common /timesteps/ dt_sub, dnt, dtmix, dtbio, dtshap

      ibio = ntrac*npt*nz + 1
      ibio1 = ntrac*npt*(nz+1) + 1
      if (use_bio_old) then
         call BIO_LIMIT_old(nstep,TR(ibio),ntrac_bio,NZ,NZI,NPT)
         call OPTIC_old(TR(ibio),ntrac_bio,H,NZ,NZI,NPT,PDENS,qsr,
     *        XZM, XZE, XPAR, XPARZE, XPARZM,zparr, zparg, zpar100) 
         call BIOMASSE_old(TR(ibio),FTR(ibio),ntrac_bio,NFLBIO,H,
     *        NZ,NZI,NPT,T,KBIO,XZM,XZE,XPAR,
     *        XPARZE,XPARZM,FLXBIO,XLE,XLM,XLT,XLN,flor)
      else
         if (i_propag_step.eq.0) then
            call DAYLENGTH (NPT, tmonth, y_deg, DAYL)
            call PROPAG(TR(ibio),ntrac_bio,H,NZ,NZI,NPT,FPAR)
         endif
               
c      in the above we are using flxbio as a work arrays
         call PAR (NPT,NZ,KBIO,tmonth,DAYL,qsr,FPAR,dtbio,XPAR)
         call ZEZM (PDENS,FPAR,H,NZ,NZI,NPT,XZM,XZE)
         call BIOMASSE(TR(ibio),FTR(ibio),ntrac_bio,NFLBIO,H,
     *        NZ,NZI,NPT,T,KBIO,XPAR,FLXBIO,XLE,XLT,XLN,flor)
         call BIO_LIMIT(nstep,TR(ibio),ntrac_bio,NZ,NZI,NPT,
     *        clm_coef_trb,trclim(ibio),inittr(ntrac+1),h,flor,dt_sub,ftr_sp(ibio1))
      endif
      i_propag_step = mod(i_propag_step+ 1, nstep_propag)

      return
      end