/***************************************************************************
                 module classe sbcgyre_qsr_qns.h  -  description
***************************************************************************/
// Mohamed Berrada
// locean-ipsl.upmc, Paris, Mars 16, 2009
// (jul) Qns : non solar heat flux (non penetrative)
// (jul) Qrs : solar heat flux (penetrative)
//===========================================================================
//                          methode forward
forward (YREAL x1)
{/*	1	from	ta    	  		1  i 	j    1    t-2*/

  /* verification code au 1er pas de temps, Mohamed, April 3, 2009
     norm(qsrY-qsrF)=     9.774031241559886e-13
     norm(qsrY-qsrF)/norm(qsrF)=     3.299060245612050e-16
     norm(qnsY-qnsF)= 4.640102485953581e-12
     norm(qnsY-qnsF)/norm(qsrF)=3.428336393249552e-16
  */
  if(Yt==TU){
    YS1=0.;
    YS2=0.;
  }
  else{
    //      ! current day (in hours) since january the 1st of the current year
    /*    double  ztime = double(Yt-TU) * rdt / (rmmss * rhhmm)  //       !  total incrementation (in hours)
	  - (nyear  - 1) * rjjhh * raajj;  //           !  minus years since beginning of experiment (in hours)
    */
    double  ztime = kt * rdt / (rmmss * rhhmm)  //       !  total incrementation (in hours)
      - (nyear  - 1) * rjjhh * raajj;  //           !  minus years since beginning of experiment (in hours)
    double  ztimemax1 = ((5.*30.)+21.)* 24.;//                     ! 21th june     at 24h in hours
    double  ztimemin1 = ztimemax1 + rjjhh * raajj / 2.;//            ! 21th december        in hours
    double  ztimemax2 = ((6.*30.)+21.)* 24.;//                      ! 21th july     at 24h in hours
    double  ztimemin2 = ztimemax2 - rjjhh * raajj / 2;//            ! 21th january         in hours
    // !                                                    ! NB: rjjhh * raajj / 4 = one seasonal cycle in hours
    
    //  ! amplitudes
    double  zTstar    = 28.3;//      ! intemsity from 28.3 a -5 deg
    
    //      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
    double  zcos_sais1 = cos( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi );
    double  zcos_sais2 = cos( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi );
    
    double  ztrp= - 40.e0;       // ! retroaction term on heat fluxes (W/m2/K)
    //  double  zconv = 3.16e-5;      //! convertion factor: 1 m/yr => 3.16e-5 mm/s
    double  t_star = zTstar*(1+1./50.*zcos_sais2)*cos(rpi*(gphit(Yi,Yj)-5.)/(53.5*(1.+11./53.5*zcos_sais2)*2.));
    //   ! 23.5 deg : tropics
    YS1 = 230*cos(3.1415*(gphit(Yi,Yj)-23.5*zcos_sais1)/(0.9*180.));
    YS2 = ztrp*(x1-t_star)-YS1; // on s'est permi ceci (YS1) car YS1 ne retropropage pas de gradient
    if(Yi==9 && Yj==1 && Yt==3){
      /*     printf("gphit(Yi,Yj)=%25.16e\n",gphit(Yi,Yj));
      printf("zcos_sais1=%25.16e\n",zcos_sais1);
      printf("YS1=%25.16e\n",YS1);*/
    }
  }
  //
}

//===========================================================================
//                         methode  backward

backward (YREAL x1)
{
  if(Yt==TU){
    YJ1I1=0.;
    YJ2I1=0.;
  }
  else{
    double  ztrp= - 40.e0;       // ! retroaction term on heat fluxes (W/m2/K)
    YJ1I1=0.;
    YJ2I1=ztrp;
  }
  //
  
}

  //===========================================================================
//********************* FIN DU MODULE sbc_gyre_qsr_qns **********************