Changes in / [20:30]


Ignore:
Location:
/trunk
Files:
8 added
6 deleted
5 edited

Legend:

Unmodified
Added
Removed

  • /trunk/condmag_from_orca.pro

    r20 r30  
    7777; 
    7878; @history 
     79; reee522 2007-06-06T15:14:24Z rhodes (IRIX64) 
     80; replace call of initorca2_bab by initocemesh 
     81; reee522 2007-06-06T14:57:39Z rhodes (IRIX64) 
     82; correction for DRAKKAR_EXP usage 
    7983; fplod 2007-03-21T13:23:55Z aedon.locean-ipsl.upmc.fr (Darwin) 
    8084; create from condmag_on_orca to validate step1 of GEOMAG 
    … …  
    138142                    ENDCASE 
    139143                    filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc' 
    140                  END 
     144                 ENDIF ELSE BEGIN 
     145                    msg = 'eee : unset DRAKKAR_EXP keyword' 
     146                    ras = report(msg) 
     147                    msg = 'eee : orcares must be G42 or G70' 
     148                    ras = report(msg) 
     149                    RETURN 
     150                 ENDELSE 
    141151              END 
    142152      ELSE  : BEGIN 
    … …  
    324334;---- 
    325335; 
    326 CASE orcares OF 
    327   'ORCA2': BEGIN 
    328               @initorca2_bab 
    329            END 
    330   ENDCASE 
     336initocemesh, orcares, DRAKKAR_EXP = drakkar_exp 
    331337; 
    332338;---- 
    … …  
    425431; ++ tri par ordre alpha , par ordre de pourcentage croissant ou décroissant 
    426432; +++ d'utilisation 
    427   profiler,/REPORT  
    428 ;  shut down all profiling (according to  
     433  profiler,/REPORT 
     434;  shut down all profiling (according to 
    429435; http://www.dfanning.com/code_tips/whyslow.html) 
    430   profiler, /CLEAR, /SYSTEM  
     436  profiler, /CLEAR, /SYSTEM 
    431437  profiler, /CLEAR, /RESET 
    432438ENDIF 

  • /trunk/condmag_on_orca.pro

    r20 r30  
    5757; 
    5858; @history 
     59; reee522 2007-06-06T14:57:39Z rhodes (IRIX64) 
     60; correction for DRAKKAR_EXP usage 
    5961; reee522 2006-12-13T13:50:32Z rhodes (IRIX64) 
    6062; add optional keyword drakkar_exp 
    … …  
    135137                    ENDCASE 
    136138                    filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc' 
    137                  END 
     139                 ENDIF ELSE BEGIN 
     140                    msg = 'eee : unset DRAKKAR_EXP keyword'  
     141                    ras = report(msg) 
     142                    msg = 'eee : orcares must be G42 or G70' 
     143                    ras = report(msg) 
     144                    RETURN 
     145                 ENDELSE 
    138146              END 
    139147      ELSE  : BEGIN 

  • /trunk/divfred.pro

    r20 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

    r20 r30  
    1 PRO forcagequimarche,iyear,ian  
    2 @init2 
    3 @initorca2_bab 
    4  
    5  
    6  
     1;+ 
     2; 
     3; @file_comments 
     4; calcul d'une matrice 3D de conductivité sigma=f(T,S) en fonction de T et S 
     5; de ORCA, sur la grille T 
     6; 
     7; @categories 
     8; orca grid 
     9; 
     10; @param orcares {in}{required}{type=string} 
     11; code of ORCA grid to use for output results 
     12; must be 'ORCA2' 
     13; 
     14; @param iyear {in}{required}{type=integer} 
     15; real year 
     16; 
     17; @param ian {in}{required}{type=integer} 
     18; simulation year 
     19; 
     20; @keyword DRAKKAR_EXP {type=string} 
     21; code for Drakkar experiment 
     22; only used when orcares = ORCA025 
     23; must be G42 ++ G70 
     24; 
     25; @restrictions 
     26;  - Requires SAXO tools 
     27; 
     28; @todo 
     29; validation of ORCA2 processing (now some nan in output files !) 
     30; add ORCA025 
     31; provide tools to plot output files 
     32; produce a NetCDF GDT or CF compliant 
     33; add divchoice parameter to choose either div or divfred 
     34; 
     35; @pre 
     36; see geomag_profile.sh 
     37; be sure to have ++ in the directory defined in ${GEOMAG_ID}/ 
     38; 
     39; @post 
     40; see geomag_profile.sh 
     41; ++ in the directory defined in ${GEOMAG_OD}/ 
     42; 
     43; @examples 
     44; for ORCA2 
     45; IDL> .run forcagequimarche 
     46; IDL> forcagequimarche, 'ORCA2', '1993', '01' 
     47; 
     48; @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) 
     51; fplod 2007-07-27T10:29:18Z cerbere.locean-ipsl.upmc.fr (Linux) 
     52; add orcares in the name of output files 
     53; fplod 2007-07-27T10:29:18Z cerbere.locean-ipsl.upmc.fr (Linux) 
     54; add orcares parameter 
     55; add optional keyword drakkar_exp 
     56; add header for idldoc 
     57; add check of io directories 
     58; replace io directories user dependant to $GEOMAG_ID and $GEOMAG_OD 
     59; replace call of divfred by call of div 
     60; 
     61; @version 
     62; $Id$ 
     63; 
     64;- 
     65; 
     66PRO forcagequimarche,orcares,iyear,ian,DRAKKAR_EXP = drakkar_exp 
     67; 
     68;---- 
     69; check input parameters 
     70;---- 
     71; 
     72; check orcares definition 
     73; 
     74  CASE orcares OF 
     75     'ORCA2': BEGIN 
     76                 msg = 'iii : valid orcares parameter = ' + orcares 
     77                 ras = report(msg) 
     78                 filename_oce = 'meshmask_bab.nc' 
     79                 IF keyword_set(DRAKKAR_EXP) THEN BEGIN 
     80                    msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp 
     81                    ras = report(msg) 
     82                 END 
     83              END 
     84      ELSE  : BEGIN 
     85                 msg = 'eee : invalid orcares parameter = ' + orcares 
     86                 ras = report(msg) 
     87                 msg = 'eee : orcares must be ORCA2 or ORCA025++' 
     88                 ras = report(msg) 
     89                 RETURN 
     90              END 
     91  ENDCASE 
     92 
     93; ++ check iyear 
     94; ++ check ian 
     95 
     96;++@init2 
     97; init grid 
     98initocemeshmask, orcares, DRAKKAR_EXP = drakkar_exp 
    799 
    8100@common 
    9  
    10  
    11  
    12 ;ian='01' 
    13 ;iyear='1993' 
    14  
     101; 
     102; check for input files 
     103; 
     104; test if ${GEOMAG_ID} defined 
     105  geomag_id_env=GETENV('GEOMAG_ID') 
     106  CASE geomag_id_env OF 
     107     ''  :  BEGIN 
     108              msg = 'eee : ${GEOMAG_ID} is not defined' 
     109              ras = report(msg) 
     110              RETURN 
     111            END 
     112     ELSE: BEGIN 
     113             msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env 
     114             ras = report(msg) 
     115           END 
     116  ENDCASE 
     117; 
     118  iodirin = isadirectory(geomag_id_env) 
     119; 
     120; existence and protection of ${GEOMAG_ID} 
     121  IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN 
     122     msg = 'eee : the directory' + iodirin  + ' is not accessible.' 
     123     ras = report(msg) 
     124     RETURN 
     125  ENDIF 
     126; 
     127; test if ${GEOMAG_OD} defined 
     128  geomag_od_env=GETENV('GEOMAG_OD') 
     129  CASE geomag_od_env OF 
     130     '' : BEGIN 
     131             msg = 'eee : ${GEOMAG_OD} is not defined' 
     132             ras = report(msg) 
     133             RETURN 
     134          END 
     135     ELSE: BEGIN 
     136             msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env 
     137             ras = report(msg) 
     138           END 
     139  ENDCASE 
     140; 
     141; check if output data will be possible 
     142  iodirout = isadirectory(geomag_od_env) 
     143; 
     144; existence and protection 
     145  IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN 
     146     msg = 'eee : the directory' + iodirout  + ' was not found.' 
     147     ras = report(msg) 
     148     RETURN 
     149  ENDIF 
     150; 
    15151e_exp='ESS' 
    16 rep_fred='/usr/work/sur/fvi/OPA/geomag/' 
    17152key_portrait = 0 
    18 ; stockage des fichiers brut  
    19   ioDATA='/usr/work/sur/fvi/OPA/ORCA2/' 
    20   file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc'  
    21   print, file_U 
    22   file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc'  
    23   print, file_V 
    24   file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc'  
    25   print, file_T 
    26   file_Sed= rep_fred+'cond_sed_ORCA2.nc' 
    27   file_Br= rep_fred+'Br_ORCA2.nc' 
    28  
    29  
     153; stockage des fichiers brut 
     154  ioDATA=geomag_id_env 
     155 
     156  CASE orcares OF 
     157     'ORCA2': BEGIN 
     158  file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 
     159  file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 
     160  file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc' 
     161  file_Sed= geomag_id_env+'cond_sed_ORCA2.nc' 
     162  file_Br= geomag_id_env+'Br_ORCA2.nc' 
     163     END 
     164  ENDCASE 
     165; 
    30166; title 
    31167     t_exp= e_exp 
    32 ;     t_exp0= e_exp0 
    33168     t_bt  = 'bar_transp' 
    34 ioORLN2 = '/usr/work/sur/fvi/OPA/ORCA2' 
     169ioORLN2 = geomag_id_env 
     170; 
    35171;facteur d'echelle vertical  for partial steps 
    36172e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 
    … …  
    39175;BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br) 
    40176BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill) 
    41  
    42  
    43  
    44  
     177; 
    45178; vertical integration: 
    46 e3t3d=make_array(jpi,jpj,jpk)    
     179e3t3d=make_array(jpi,jpj,jpk) 
    47180        for k=0, jpk-1 do begin              &$ 
    48           for j=0,jpj-1 do begin              &$       
     181          for j=0,jpj-1 do begin              &$ 
    49182            for i=0,jpi-1 do begin             &$ 
    50183              e3t3d(i,j,k) = e3t(k)    &$ 
    … …  
    53186        endfor 
    54187jpt = 73 
    55  
    56 ;vud = make_array(jpi,jpj,jpt)    
     188; 
     189;vud = make_array(jpi,jpj,jpt) 
    57190;vvd = make_array(jpi, jpj, jpt) 
    58191divBustar = make_array(jpi, jpj, jpt) 
    59192diver3=replicate(0.,182,149,1,73) 
    60193sigma3=replicate(0.,182,149,1,73) 
    61  
    62  
    63 FOR jt = 0, jpt-1 DO BEGIN &$    
    64  
     194; 
     195FOR jt = 0, jpt-1 DO BEGIN &$ 
     196; 
    65197; ouverture des fichiers et stockage en memoire partial steps 
    66   vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U)  &$    
    67 ;stop 
    68   vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V)  &$    
     198  vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U)  &$ 
     199;stop 
     200  vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V)  &$ 
    69201;stop 
    70202; lecture salinite & temperature 
    … …  
    77209  conduct(mask_t)=0. 
    78210; Somme conduct au point T 
    79  
    80  
     211; 
    81212; Calcul SIGMA 
    82  
     213; 
    83214  SIGMAoc=total(conduct*e3t3d*tmask,3) 
    84215  SIGMA=SIGMAsed+SIGMAoc 
    85  
    86  
    87  
     216; 
    88217  SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2. 
    89  
    90218  SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2. 
    91  
     219; 
    92220; Calcul B en points u et v 
    93  
     221; 
    94222  BR_u=(BR+shift(BR,-1,0))/2. 
    95  
    96223  BR_v=(BR+shift(BR,0,1))/2. 
    97  
    98  
     224; 
    99225; Calcul integrale conduct 
    100  
     226; 
    101227  conduct_u=(conduct+shift(conduct,-1,0,0))/2. 
    102  
    103228  conduct_v=(conduct+shift(conduct,0,1,0))/2. 
    104    
     229; 
    105230  u_cond_u=total( vu*conduct_u*e3u3d*umask(),3) 
    106  
    107231  v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3) 
    108  
    109  
     232; 
    110233  Bu_star= BR_u*u_cond_u/SIGMA_u 
    111  
    112234  Bv_star= BR_v*v_cond_v/SIGMA_v 
    113  
     235; 
    114236; Divergence du champ 
    115  
    116237  Diver=divfred(Bu_star,Bv_star) 
    117238 
    … …  
    120241lecontinent=where(Diver GT 1.E08) 
    121242Diver(lecontinent)=0. 
    122  
    123  
    124243 
    125244;stop 
    … …  
    134253sigma3(181,*,*,*)=sigma3(1,*,*,*) 
    135254 
    136  
    137255print,  jt 
    138256 
    139  
    140 ENDFOR   
     257ENDFOR 
    141258; on ferme le NetCDF 
    142259;NCDF_CLOSE,id3 
    143260;NCDF_CLOSE,id4 
    144  
    145  
    146  
    147261 
    148262temps=fltarr(73) 
    … …  
    154268 
    155269   vargrid = 'T' 
    156    iodir = '/usr/work/sur/fvi/OPA/ORCA2/' 
     270   iodir = geomag_od_env 
    157271; Nom 
    158    idout = NCDF_CREATE(iodir+'DivBustar_5d_'+iyear+'_grid_T.nc',/clobber) 
     272   idout = NCDF_CREATE(iodir+'DivBustar_5d_'+iyear+'_grid_T_'+orcares+'.nc',/clobber) 
    159273   print, 'Creation du fichier Netcdf' 
    160274   NCDF_CONTROL, idout, /nofill 
    … …  
    203317; Creation de la latitude 
    204318   NCDF_VARPUT, idout, id1, gphit 
    205 ; Creation de la profondeur  
    206    NCDF_VARPUT, idout, id2, gdept  
     319; Creation de la profondeur 
     320   NCDF_VARPUT, idout, id2, gdept 
    207321; Creation du calendrier 
    208322 
    209323   NCDF_VARPUT, idout, id3, temps 
    210324 
    211     
     325 
    212326; Ecriture des donnees 
    213327 
    214 ; ecriture des glam_8  
     328; ecriture des glam_8 
    215329   NCDF_VARPUT, idout, id4 , diver3 
    216330 
    … …  
    220334;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
    221335 
    222   idout = NCDF_CREATE(iodir+'Sigma_5d_'+iyear+'_grid_T.nc',/clobber) 
     336  idout = NCDF_CREATE(iodir+'Sigma_5d_'+iyear+'_grid_T_'+orcares+'.nc',/clobber) 
    223337   print, 'Creation du fichier Netcdf' 
    224338   NCDF_CONTROL, idout, /nofill 
    … …  
    267381; Creation de la latitude 
    268382   NCDF_VARPUT, idout, id1, gphit 
    269 ; Creation de la profondeur  
    270    NCDF_VARPUT, idout, id2, gdept  
     383; Creation de la profondeur 
     384   NCDF_VARPUT, idout, id2, gdept 
    271385; Creation du calendrier 
    272386 
    273387   NCDF_VARPUT, idout, id3, temps 
    274388 
    275     
     389 
    276390; Ecriture des donnees 
    277391 
    278 ; ecriture des glam_8  
     392; ecriture des glam_8 
    279393   NCDF_VARPUT, idout, id4 , sigma3 
    280394 

  • /trunk/geomag_profile.sh

    r20 r30  
    66# usage : 
    77# online usage 
    8 # $ . geomag_profile.sh -d directory 
     8# $ . geomag_profile.sh -d directory -i indir -o outdir -t tempdir 
    99# 
    1010# in ${HOME}/.profile, add the following line 
    11 # . geomag_profile.sh -d directory 
     11# . geomag_profile.sh -d directory  -i indir -o outdir -t tempdir 
    1212# 
    1313# examples : 
    14 # for fplod, on aedon.locean-ipsl.upmc.fr 
     14# for fplod, on aedon.locean-ipsl.upmc.fr  
    1515# $ cd /usr/home/fplod/incas/geomag/geomag_ws 
    16 # $  . geomag_profile.sh -d /usr/home/fplod/incas/geomag/geomag_ws/ 
     16# $  . geomag_profile.sh -d /usr/home/fplod/incas/geomag/geomag_ws/ -i ${HOME}/geomag_d/ -o ${HOME}/geomag_d/ -t /usr/temp/${LOGNAME}/log/ 
     17# 
    1718# for reee522 on rhodes.idris.fr 
    1819# $ cd ${HOME}/incas/geomag/geomag_ws/ 
    19 # $ . geomag_profile.sh -d ${HOME}/incas/geomag/geomag_ws/ 
     20# $ . geomag_profile.sh -d ${HOME}/incas/geomag/geomag_ws/ -i ${HOMEGAYA}/geomag_d/ -o ${HOMEGAYA}/geomag_d/ -t /tmp/${LOGNAME}/log/ 
    2021# 
    2122# original location 
    … …  
    2627# ++ machine dependant 
    2728# ++ besoin de posix 
    28 # ++ pas de MANPATH defini par déut sur zeus 
     29# ++ pas de MANPATH defini par déubt sur zeus 
     30# reee522 2007-06-14T14:25:00Z rhodes (IRIX64) 
     31# add -i, -o and -t parameters to be more flexible 
     32# some bug fixes 
     33# more checking 
    2934# fplod 2007-03-13T16:12:41Z zeus.locean-ipsl.upmc.fr (Linux) 
    3035# add creation of $GEOMAG_ID and $GEOMAG_OD if they don't exist 
    … …  
    3843command=geomag_profile.sh 
    3944# 
    40 usage=" Usage : ${command} -d directory" 
     45usage=" Usage : ${command} -d directory -i indir -o outdir -t tempdir" 
    4146# 
    4247while [ ! -z "${1}" ] # ++ pb bash 
    4348do 
    4449 case ${1} in 
    45  -d) # directory choosen by user 
     50 -d) # directory for application choosen by user (see svn checkout command used) 
    4651  directory=${2} 
    4752  shift 
    4853 ;; 
     54 -i) # directory for inputs choosen by user 
     55  indir=${2} 
     56  shift 
     57 ;; 
     58 -o) # directory for outputs choosen by user 
     59  outdir=${2} 
     60  shift 
     61 ;; 
     62 -t) # directory for temporary outputs choosen by user 
     63  tempdir=${2} 
     64  shift 
     65 ;; 
     66 
    4967 *) # other choice 
    5068  echo "${usage}" 
    … …  
    7189set -u 
    7290system=$(uname) 
    73 case ${system} in 
     91case "${system}" in 
    7492IRIX64) 
    7593  echo " www : no specific posix checking" 
    … …  
    7997;; 
    8098esac 
    81 # 
    8299# 
    83100GEOMAG=${directory} 
    … …  
    110127fi 
    111128# 
    112 GEOMAG_LOG=/usr/temp/${LOGNAME}/log/ 
     129GEOMAG_LOG=${tempdir} 
    113130export GEOMAG_LOG 
     131if [ ! -d ${GEOMAG_LOG} ]  
     132then 
     133  mkdir -p ${GEOMAG_LOG} 
     134  echo "${command} : iii : creation of \${GEOMAG_LOG}" 
     135fi  
     136# check for permission on GEOMAG_LOG 
     137if [ ! -x ${GEOMAG_LOG} ] 
     138then 
     139 echo " eee : ${GEOMAG_LOG} not reachable" 
     140 # nb : no exit because this file should be launched by login process 
     141fi 
     142# 
     143# check for permission on GEOMAG_LOG 
     144if [ ! -w ${GEOMAG_LOG} ] 
     145then 
     146 echo " eee : ${GEOMAG_LOG} not writable" 
     147 # nb : no exit because this file shouldreachable be launched by login process 
     148fi 
    114149# 
    115150EDITOR=vi 
    … …  
    117152# 
    118153# io directories 
    119 GEOMAG_ID=${HOME}/geomag_d/ 
     154GEOMAG_ID=${indir} 
    120155export GEOMAG_ID 
    121156if [ ! -d ${GEOMAG_ID} ]  
    … …  
    124159  echo "${command} : iii : creation of \${GEOMAG_ID}" 
    125160fi  
    126 GEOMAG_OD=${HOME}/geomag_d/ 
     161# check for permission on GEOMAG_ID 
     162if [ ! -x ${GEOMAG_ID} ] 
     163then 
     164 echo " eee : ${GEOMAG_ID} not reachable" 
     165 # nb : no exit because this file should be launched by login process 
     166fi 
     167# 
     168GEOMAG_OD=${outdir} 
    127169export GEOMAG_OD 
    128 if [ ! -d ${GEOMAG_ID} ] 
     170if [ ! -d ${GEOMAG_OD} ] 
    129171then 
    130172  mkdir -p ${GEOMAG_OD} 
    131173  echo "${command} : iii : creation of \${GEOMAG_OD}" 
    132174fi 
     175# check for permission on GEOMAG_OD 
     176if [ ! -x ${GEOMAG_OD} ] 
     177then 
     178 echo " eee : ${GEOMAG_OD} not reachable" 
     179 # nb : no exit because this file should be launched by login process 
     180fi 
     181if [ ! -w ${GEOMAG_OD} ] 
     182then 
     183 echo " eee : ${GEOMAG_OD} not writable" 
     184 # nb : no exit because this file should be launched by login process 
     185fi 
    133186# 
    134187# end 
Note: See TracChangeset for help on using the changeset viewer.