;**********************************************************
; Computes and removes the seasonal cycle of a variable
; var     : 1,2,3 or 4D field
; mes     : month vector (1:12 x nyears)
; n       : dimension corresponding to time
; varm    : seasonal cycle
; var_anom: var anomalies (after removing seasonal cycle)
;**********************************************************


pro remove_seas,var,mes,n,varm,var_anom


dim = size(var)

;======================================================
; 1D VARIABLE
;======================================================

if n eq 1 then begin
	varm = fltarr(12)
	
	for imes = 1,12 do begin
		in = where(mes eq imes)
        varm(imes-1) = (total(var(in),1))/n_elements(in)
   	endfor

	var_anom = fltarr(dim(1))
   	for it = 0,dim(1)-1 do begin
       	imes = mes(it)
       	var_anom(it) = var(it) - varm(imes-1)
   	endfor
endif


;======================================================
; 2D VARIABLE
;======================================================

if n eq 2 then begin
	varm = fltarr(dim(1),12)
   	for imes = 1,12 do begin
      	in = where(month eq imes)
        varm(*,imes-1) = (total(var(*,in),2))/n_elements(in)
   	endfor

	var_anom = fltarr(dim(1),dim(2))
   	for it = 0,dim(2)-1 do begin
      	imes = mes(it)
        var_anom(*,it) = var(*,it) - varm(*,imes-1)
   	endfor
endif


;======================================================
; 3D VARIABLE
;======================================================

if n eq 3 then begin
	varm = fltarr(dim(1),dim(2),12)
   	for imes = 1,12 do begin
		in = where(mes eq imes);
        varm(*,*,imes-1) = (total(var(*,*,in),3))/n_elements(in)
   	endfor

	var_anom = fltarr(dim(1),dim(2),dim(3))
   	for it = 0,dim(3)-1 do begin
   		imes = mes(it)
        var_anom(*,*,it) = var(*,*,it) - varm(*,*,imes-1)
   	endfor
endif


;======================================================
; 4D VARIABLE
;======================================================

if n eq 4 then begin
	varm = fltarr(dim(1),dim(2),dim(3),12)
   	for imes = 1,12 do begin
      	in = where(mes eq imes);
      	varm(*,*,*,imes-1) = (total(var(*,*,*,in),4))/n_elements(in)
    endfor

	var_anom = fltarr(dim(1),dim(2),dim(3),dim(4))
   	for it = 0,dim(4)-1 do begin
   		imes = mes(it)
        var_anom(*,*,*,it) = var(*,*,*,it) - varm(*,*,*,imes-1)
   	endfor
endif

end