PRO forcage,iyear,ian 
@init2
@initorca2_bab




@common



ian='01'
iyear='1993'

e_exp='ESS'
rep_fred='/usr/work/sur/fvi/OPA/geomag/'
key_portrait = 0
; stockage des fichiers brut 
  ioDATA='/usr/work/sur/fvi/OPA/ORCA2/'
  file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 
  print, file_U
  file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 
  print, file_V
  file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc' 
  print, file_T
  file_Sed= rep_fred+'cond_sed_ORCA2.nc'
  file_Br= rep_fred+'Br_ORCA2.nc'


; title
     t_exp= e_exp
;     t_exp0= e_exp0
     t_bt  = 'bar_transp'
ioORLN2 = '/usr/work/sur/fvi/OPA/ORCA2'
;facteur d'echelle vertical  for partial steps
e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
SIGMAsed=read_ncdf('cond_sed',0,/timestep,/nostruct,/tout,filename=file_Sed,/cont_nofill)
;BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br)
BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill)




; vertical integration:
e3t3d=make_array(jpi,jpj,jpk)   
        for k=0, jpk-1 do begin              &$
          for j=0,jpj-1 do begin              &$      
            for i=0,jpi-1 do begin             &$
              e3t3d(i,j,k) = e3t(k)    &$
            endfor                     &$
          endfor                      &$
        endfor
jpt = 73

;vud = make_array(jpi,jpj,jpt)   
;vvd = make_array(jpi, jpj, jpt)
divBustar = make_array(jpi, jpj, jpt)

; ouverture des fichiers dans lesquels on va écrire 
;id3=NCDF_OPEN('/usr/work/sur/fvi/OPA/geomag/U_5d_'+iyear+'_grid_T.nc',/write)
;id4=NCDF_OPEN('/usr/work/sur/fvi/OPA/geomag/V_5d_'+iyear+'_grid_T.nc',/write)
id4=NCDF_CREATE('/usr/work/sur/fvi/OPA/ORCA2/DivBustar_5d_'+iyear+'_grid_T.nc',/clobber)

;id4=NCDF_OPEN('/usr/work/sur/fvi/OPA/ORCA2/DivBustar_5d_'+iyear+'_grid_T.nc',/write)

FOR jt = 0, jpt-1 DO BEGIN &$   

; ouverture des fichiers et stockage en memoire partial steps
  vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U)  &$   
  vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V)  &$   
; lecture salinite & temperature
  temp= read_ncdf('votemper',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T)
  salin=read_ncdf('vosaline',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T)

  conduct=0.02047780622061 + 0.00273147624197*temp + 0.00035133182334*temp*temp + 0.09139808809909*salin + 0.00241425798890*salin*temp -0.00023998958774*salin*salin
  mask_t=where(conduct GT 1.e+19)
  conduct(mask_t)=0.
; Somme conduct au point T


; Calcul SIGMA

  SIGMAoc=total(conduct*e3t3d*tmask,3)
  SIGMA=SIGMAsed+SIGMAoc

  SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2.

  SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2.

; Calcul B en points u et v

  BR_u=(BR+shift(BR,-1,0))/2.

  BR_v=(BR+shift(BR,0,1))/2.


; Calcul integrale conduct

  conduct_u=(conduct+shift(conduct,-1,0,0))/2.

  conduct_v=(conduct+shift(conduct,0,1,0))/2.
  
  u_cond_u=total( vu*conduct_u*e3u3d*umask(),3)

  v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3)


  Bu_star= BR_u*u_cond_u/SIGMA_u

  Bv_star= BR_v*v_cond_v/SIGMA_v

; Divergence du champ

  Diver=div(Bu_star,Bv_star)
  
; Somme sur la verticale partial steps
;                vum=total( vu*e3u3d*umask(),3 )  &$   
;                vvm=total( vv*e3v3d*vmask(),3 )  &$   

; Shift sur la grille T partial steps
;             vut= (vum+shift(vum,1,0) )*0.5   &$   
;             vvt= (vvm+shift(vvm,0,1) )*0.5   &$
; Bande de recouvrement 
;             vut(0, *) = vut(jpi-2, *)
;             vvt(*, 0) = 0.
; stockage dans le fichier de sortie 
;NCDF_VARPUT, id3,'sossheig',vut, offset = [0, 0, jt]
NCDF_VARPUT, id4,'sossheig',Diver, offset = [0, 0, jt]

print,  jt

ENDFOR  
; on ferme le NetCDF
;NCDF_CLOSE,id3
NCDF_CLOSE,id4




END