Changeset 30


Ignore:
Timestamp:
07/30/07 15:27:55 (17 years ago)
Author:
pinsard
Message:

add step2_diff.pro (not ok yes) and restore divfred call

Location:
trunk
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/divfred.pro

    r12 r30  
    1 ;------------------------------------------------------------ 
    2 ;------------------------------------------------------------ 
    3 ;------------------------------------------------------------ 
    41;+ 
    5 ; NAME:div 
    6 ; 
    7 ; PURPOSE:calcule la divergence d'un champ 2D 
    8 ; 
    9 ; CATEGORY:calcule sur les matrices 
    10 ; 
    11 ; CALLING SEQUENCE:res=div(u,v) 
    12 ; 
    13 ; INPUTS: 
    14 ;       u et v deux matrices representant les coordonnes d''un 
    15 ;       champ de vecteur 
    16 ; 
    17 ; KEYWORD PARAMETERS: 
    18 ; 
    19 ; OUTPUTS:res: une matrice 2d 
    20 ; 
    21 ; COMMON BLOCKS: 
    22 ;       common.pro 
    23 ; 
    24 ; SIDE EFFECTS: 
    25 ; 
    26 ; RESTRICTIONS: 
    27 ; les matrices u et v peuvent de 2 a 4 dimensions. 
    28 ; attention pour distinger les differents configurations de u et v 
    29 ; (xy, xyz, xyt, xyzt), on regarde la variable du common  
     2; @files_comments 
     3; calcule la divergence d'un champ 2D 
     4; 
     5; @categories 
     6; matrix 
     7; 
     8; @examples 
     9; IDL> res=Divfred(u,v) 
     10; 
     11; @param u {in}{required}{type=2D array} 
     12; Matrix representing the zonal coordinates (U point) 
     13; 
     14; @param v {in}{required}{type=2D array} 
     15; Matrix representing the zonal coordinates (V point) 
     16; 
     17; @returns 
     18; the divergence of the input data (with the same size) 
     19; 
     20; @uses 
     21; common.pro 
     22; 
     23; @restrictions 
     24; - Requires SAXO tools 
     25; - les matrices u et v peuvent de 2 a 4 dimensions. 
     26; attention pour distinguer les différentes configurations de u et v 
     27; (xy, xyz, xyt, xyzt), on regarde la variable du common 
    3028;        -time qui contient le calendrier en jour julien d''IDL auquel 
    31 ;        se rapportent u et v ansi que la variable  
    32 ;        -jpt qui est le nombre de pas de temps a considerer ds time. 
    33 ; les tableaux u et v sont decoupes sur le meme domaine 
    34 ; geographique. A cause du decalage des grilles T, U, V et F il est 
    35 ; possiible que ces 2 tableaux n''aient pas la meme taille et se 
    36 ; repportent a des indices differents. Si tel est le cas les tableaux 
    37 ; sont redecoupes sur les indices qu'ils ont en commun et le dommaine 
    38 ; est redefinit pour qu'il colle a ces indices communs. 
    39 ; pour eviter ces redecoupes utiliser le mot cles /memeindice ds 
    40 ; domdef.pro 
    41 ; 
    42 ; les points sur le bord du dessin sont mis a !values.f_nan 
    43 ; 
    44 ; EXAMPLE: 
    45 ; 
    46 ; MODIFICATION HISTORY:Guillaume Roullet (grlod\@ipsl.jussieu.fr) 
     29;        se rapportent u et v ansi que la variable 
     30;        -jpt qui est le nombre de pas de temps à considérer dans time. 
     31; les tableaux u et v sont découpés sur le même domaine 
     32; géographique. A cause du décalage des grilles T, U, V et F il est 
     33; possible que ces 2 tableaux n''aient pas la même taille et se 
     34; reportent à des indices différents. Si tel est le cas les tableaux 
     35; sont redécoupés sur les indices qu'ils ont en commun et le dommaine 
     36; est redéfini pour qu'il colle à ces indices communs. 
     37; pour éviter ces redécoupes utiliser le mot clé /memeindice dans 
     38; <pro>domdef</pro> 
     39; 
     40; les points sur le bord du dessin sont mis à !values.f_nan 
     41; 
     42; @todo 
     43; ++ pas fini de comprendre, tester (compatibilité saxo), adapter, commenter 
     44; nettoyer 
     45; ++ a comparer et merger avec SAXO_DIR/ToBeReviewed/CALCULS/div.pro 
     46;       Written by:   fplod 2006-03-27T13:45:14Z aedon.lodyc.jussieu.fr (Darwin) 
     47; 
     48; @history 
     49; fplod 2007-07-30T10:11:33Z aedon.locean-ipsl.upmc.fr (Darwin) 
     50; add header for idldoc 
     51; created from 
     52; /usr/work/sur/fvi/OPA/geomag/divfred.pro 
     53; Guillaume Roullet (grlod\@ipsl.jussieu.fr) 
    4754;                      Creation : printemps 1998 
    4855;                      Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    4956;                      adaptation pour marcher avec un domaine reduit 
    5057;                      12/1/2000 
     58; @version 
     59; $Id$ 
     60; 
    5161;- 
    52 ;------------------------------------------------------------ 
    53 ;------------------------------------------------------------ 
    54 ;------------------------------------------------------------ 
     62; 
    5563FUNCTION divfred, uu, vv 
    5664   tempsun = systime(1)         ; pour key_performance 
     
    6674; on trouve les points que u et v ont en communs 
    6775;------------------------------------------------------------ 
    68    indicexu = (lindgen(jpi))[premierxu:premierxu+nxu-1] 
    69    indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1] 
     76   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
     77   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
    7078   indicex = inter(indicexu, indicexv) 
    71    indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1] 
    72    indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1] 
     79   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
     80   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
    7381   indicey = inter(indiceyu, indiceyv) 
    74    nx = n_elements(indicex)  
     82   nx = n_elements(indicex) 
    7583   ny = n_elements(indicey) 
    7684   indice2d = lindgen(jpi, jpj) 
     
    8189;xyz 
    8290;---------------------------------------------------------------------------- 
    83       (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN  
     91      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 
    8492;------------------------------------------------------------ 
    8593; extraction de u et v sur le domaine qui convient 
     
    8896            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1 
    8997            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    90              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
     98             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
    9199               case 1 of 
    92                   nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    93                   nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    94                   nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    95                   nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     100                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     101                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     102                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     103                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    96104                  ELSE : 
    97105               endcase 
    98106            END 
    99107            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    100              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
     108             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
    101109               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    102110               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    103111            END 
    104             ELSE:BEGIN  
     112            ELSE:BEGIN 
    105113               zdiv = -1 
    106114               GOTO, sortie 
    107115            end 
    108              
     116 
    109117         endcase 
    110118;------------------------------------------------------------ 
     
    114122         zu = reform(zu, nx, ny, nzt, /over) 
    115123         zu = temporary(u)*temporary(zu) $ 
    116           *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     124          *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    117125         terreu = where(zu EQ 0) 
    118126;;;         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan 
     
    121129         zv = reform(zv, nx, ny, nzt, /over) 
    122130         zv = temporary(v)*temporary(zv) $ 
    123           *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     131          *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    124132         terrev = where(zv EQ 0) 
    125133;;;         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan 
     
    129137         zdiv = reform(zdiv, nx, ny, nzt, /over) 
    130138         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $ 
    131           *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     139          *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    132140;------------------------------------------------------------ 
    133141; mise a !values.f_nan de la bordure 
    134142;------------------------------------------------------------ 
    135          if  NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     143         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    136144            zdiv(0, *, *) = !values.f_nan 
    137145            zdiv(nx-1, *, *) = !values.f_nan 
     
    143151; 
    144152         if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    145          terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] EQ 0) 
     153         terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0) 
    146154         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask 
    147155;------------------------------------------------------------ 
     
    151159         varname = 'div' 
    152160         varunits = '1e6*s-1' 
    153          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t'] 
     161         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 
    154162         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan, boite = boite) 
    155163      END 
     
    159167;---------------------------------------------------------------------------- 
    160168;---------------------------------------------------------------------------- 
    161       date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN  
     169      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
    162170;------------------------------------------------------------ 
    163171; extraction de u et v sur le domaine qui convient 
     
    166174            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1 
    167175            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    168              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
     176             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
    169177               case 1 of 
    170                   nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    171                   nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    172                   nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    173                   nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     178                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     179                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     180                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     181                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    174182                  ELSE : 
    175183               endcase 
    176184            END 
    177185            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    178              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
     186             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
    179187               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
    180188               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 
     
    185193; calcul de la divergence 
    186194;------------------------------------------------------------ 
    187          zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt] 
     195         zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    188196         terreu = where(zu EQ 0) 
    189197;;;         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan 
     
    192200         zu = temporary(u)*temporary(zu) 
    193201; 
    194          zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt] 
     202         zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    195203         terrev = where(zv EQ 0) 
    196204;;;         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan 
     
    199207         zv = temporary(v)*temporary(zv) 
    200208; 
    201          zdiv = 1e6*tmask[indice2d+jpi*jpj*premierzt]/(e1t[indice2d]*e2t[indice2d]) 
     209         zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d]) 
    202210         zdiv = (zdiv)[*]#replicate(1, jpt) 
    203211         zdiv = reform(zdiv, nx, ny, jpt, /over) 
     
    207215; mise a !values.f_nan de la bordure 
    208216;------------------------------------------------------------ 
    209          if  NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     217         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    210218            zdiv(0, *, *) = !values.f_nan 
    211219            zdiv(nx-1, *, *) = !values.f_nan 
     
    222230         varname = 'div' 
    223231         varunits = '1e6*s-1' 
    224          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t'] 
     232         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 
    225233         if keyword_set(direc) then  zdiv = grossemoyenne(zdiv,direc,/nan, boite = boite) 
    226234      END 
     
    230238;---------------------------------------------------------------------------- 
    231239;---------------------------------------------------------------------------- 
    232       date1 NE date2 AND (size(u))[0] EQ 4:BEGIN  
     240      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 
    233241         return, report('non code!') 
    234242      END 
     
    240248      ELSE:BEGIN                ;xy 
    241249         indice3d = lindgen(jpi, jpj, jpk) 
    242          indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, niveau -1] 
     250         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt] 
    243251;------------------------------------------------------------ 
    244252; extraction de u et v sur le domaine qui convient 
    245253;------------------------------------------------------------ 
    246254         case 1 of 
    247             (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN  
     255            (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN 
    248256               zdiv = -1 
    249257               GOTO, sortie 
    250258            end 
    251259            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 
    252              (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
     260             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
    253261               case 1 of 
    254                   nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
    255                   nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
    256                   nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
    257                   nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
     262                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
     263                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
     264                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
     265                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
    258266                  ELSE : 
    259267               endcase 
    260268            END 
    261269            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
    262              (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN  
     270             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
    263271               u = u[indice2d] 
    264272               v = v[indice2d] 
     
    280288; mise a !values.f_nan de la bordure 
    281289;------------------------------------------------------------ 
    282          if  NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     290         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    283291            zdiv(0, *) = !values.f_nan 
    284292            zdiv(nx-1, *) = !values.f_nan 
     
    298306         varname = 'div' 
    299307         varunits = '1e6*s-1' 
    300          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t'] 
     308         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 
    301309         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan, boite = boite) 
    302310 
    303311      END 
    304312   endcase 
    305     
     313 
    306314;------------------------------------------------------------ 
    307315sortie: 
    308    if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun  
    309     
     316   if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun 
     317 
    310318   return, zdiv 
    311319end 
  • trunk/forcagequimarche.pro

    r29 r30  
    1 ; + 
     1;+ 
    22; 
    33; @file_comments 
     
    3131; provide tools to plot output files 
    3232; produce a NetCDF GDT or CF compliant 
     33; add divchoice parameter to choose either div or divfred 
    3334; 
    3435; @pre 
     
    4647; 
    4748; @history 
     49; fplod 2007-07-30T10:12:38Z aedon.locean-ipsl.upmc.fr (Darwin) 
     50; restore divfred call instead of div ("official" saxo divergence module) 
    4851; fplod 2007-07-27T10:29:18Z cerbere.locean-ipsl.upmc.fr (Linux) 
    4952; add orcares in the name of output files 
     
    232235; 
    233236; Divergence du champ 
    234   Diver=div(Bu_star,Bv_star) 
     237  Diver=divfred(Bu_star,Bv_star) 
    235238 
    236239Diver=Diver*1e-6 
Note: See TracChangeset for help on using the changeset viewer.