source: trunk/tools/bsf/build_laplacien_f.pro @ 2

;+
;;   
;; 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
;;
;; 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
;;
;-
PRO build_laplacien_f
@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