source: trunk/tools/bsf/build_laplacien_f.pro

;+
;   
; Construction du laplacien 2D
; cherche a le lire dans le fichier
; '/workdir/rech/ces/rces023/laplacien_T.dat', sinon il le calcul   
;
; le laplacien correspond a l''ancienne grille globale d''OPA
;
; il faut modifier cette routine pour toute autre grille
;
; @uses
; <pro>common</pro>
; <propost_it>com_eg</propost_it>
;
; @history
; Creation     : juin 1998 :  G. Roullet
; Modification : 26/02/00  :  G.Roullet 
;                                Adpation pour ORCA : repliement au pole
;                                Changement du traitement des conditions
;                                Est-Ouest (on masque au lieu de supprimer les points)
;                                   => routine inversion corrigee
;
; @version
; $Id$
;
;-
PRO build_laplacien_f
;
  compile_opt idl2, strictarrsubs
;
@common
@com_eg
;
   COMMON laplaciens, lapla_t, lapla_f

   N = jpi
   N2 = long(N)*(jpj)
   bande = dblarr(N2, 5)


   fm = fmaskr[*, *, 0]/(e1f*e2f)

;
; On masque le domaine non physique (points en double)
;
   fm[0, *] = 0
   fm[jpi-1, *] = 0
;vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
; specificite ORCA
   fm[*, jpj-1] = 0
   fm[*, jpj-2] = 0
;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   um = umaskr[*, *, 0]*e1u/e2u
   vm = vmaskr[*, *, 0]*e2v/e1v
;
;  construction des 5 bandes
;

;      Sud
   bande[ 0:N2-1 , 1] = fm*shift(um, 0, 0)
;      Ouest                 
   bande[ 0:N2-1 , 2] = fm*shift(vm, 0, 0)
;      Est                   
   bande[ 0:N2-1 , 3] = fm*shift(vm,-1, 0)
;      Nord                  
   bande[ 0:N2-1 , 4] = fm*shift(um, 0,-1)
;
;
;
   bande[0:N-1      , 1] = 0.
;
;  elements ci dessous non nuls du fait des CL periodiques E-W
;  decommenter les 2 lignes si votre domaine est ferme a l''Est et a l''Ouest
;   bande[0          , 2] = 0.
;   bande[N2-1       , 3] = 0.

   bande[N2-N:N2-1  , 4] = 0.
   
;      Centre
   bande[ * , 0] = - bande[*, 1] -  bande[*, 2] -  bande[*, 3] -  bande[*, 4]

;  
;  construction de la matrice creuse
;

   ;; nb d'elements de la matrice creuse
   Nb = N2+n_elements(where(bande[*, 1:4] NE 0))+1

   ;; definition de la matrice creuse
   lapla_f = { laplacien_f, sa:dblarr(Nb), ija:lonarr(Nb) }

   ;; elements diagonaux
   lapla_f.sa[0:N2-1] = bande[0:N2-1, 0]
   lapla_f.sa[N2] = 0
   lapla_f.ija[N2] = Nb+1
   lapla_f.ija[0] = N2+2

;   
; elements extra-diagonaux
;
   ;;  -  i balaye dans l'ordre tout les elements susceptibles d'etre
   ;;     non nuls dans l'immense matrice i.e. 4 colonnes par ligne.

   ;;  -  k repere l'indice de la matrice creuse courant

   ;;  -  ligne et colonne sont les coordonnees dans la matrice pleine

   ;;  -  zligne et zcol reperent la position du point courant dans la
   ;;          matrice jpixjpj

   ;;  - decalage repere la position *relative* des 4 voisins dans la matrice
   ;;    creuse                       ^^^^^^^^

   i = long[0]
   k = long[N2]
   compteur = long[0]
   REPEAT BEGIN 
      ligne = i/4
      ;;
      ;; conditions aux limites periodiques selon W-E
      ;;
      
      zcol = ligne MOD N
      zligne = ligne/N

      IF i MOD 4 EQ 0 THEN begin
         decalage = [-N, -1, 1, N]
         ordre = [1, 2, 3, 4]
         IF zcol EQ 1 THEN BEGIN               
            decalage = [-N, -3+N, 1, N]
            ordre = [1, 3, 2, 4]
         ENDIF 
         IF zcol EQ (N-2) THEN BEGIN 
            decalage = [-N, -1, +3-N, N]
            ordre = [1, 3, 2, 4]
         ENDIF 
;vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
;
; ajout ORCA
;
; conditions de repliement au Nord
;
         IF (zligne EQ jpj-3) THEN BEGIN
            ptN = 2*(N/2-zcol)-1
            decalage[3] = ptN
            IF zcol GE N/2 THEN ordre = [1, 4, 2, 3]
         ENDIF 
;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      ENDIF 

      colonne = ligne+decalage[ordre[i MOD 4]-1]
;
      IF (colonne GE 0) AND (colonne LT N2) THEN  BEGIN 
         x = bande[ligne, ordre[i MOD 4] ]
         IF x NE 0. THEN BEGIN 
            k = k+1
            lapla_f.sa[k] = x
            lapla_f.ija[k] = colonne+1
         ENDIF ELSE BEGIN 
         ENDELSE 
      ENDIF
;
      IF (i MOD 4) EQ 3 THEN lapla_f.ija[ligne+1] = k+2
;
      i = i+1
   ENDREP UNTIL i EQ 4*N2-1

END