Ignore:
Timestamp:
05/02/06 14:59:12 (18 years ago)
Author:
pinsard
Message:

upgrade of CALCULS according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files

File:
1 copied

Legend:

Unmodified
Added
Removed
  • trunk/ToBeReviewed/CALCULS/curl.pro

    r23 r25  
    6161   tempsun = systime(1)         ; pour key_performance 
    6262; 
     63   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 
     64     return, report(['This version of curl is based on Arakawa C-grid.' $ 
     65                     , 'U and V grids must therefore be defined']) 
     66; 
    6367   u = litchamp(uu) 
    6468   v = litchamp(vv) 
     
    7175; on trouve les points que u et v ont en communs 
    7276;------------------------------------------------------------ 
    73    indicexu = (lindgen(jpi))[premierxu:premierxu+nxu-1] 
    74    indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1] 
     77   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 
     78   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 
    7579   indicex = inter(indicexu, indicexv) 
    76    indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1] 
    77    indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1] 
     80   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 
     81   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 
    7882   indicey = inter(indiceyu, indiceyv) 
    7983   nx = n_elements(indicex)  
     
    9599             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    96100            case 1 of 
    97                nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
    98                nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    99                nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
    100                nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     101               nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
     102               nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     103               nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
     104               nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    101105               ELSE : 
    102106            endcase 
     
    114118         coefu = (e1u[indice2d])[*]#replicate(1, nzt) 
    115119         coefu = reform(coefu, nx, ny, nzt, /over) 
    116          coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     120         coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    117121         terreu = where(coefu EQ 0) 
    118122         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
     
    120124         coefv = (e2v[indice2d])[*]#replicate(1, nzt) 
    121125         coefv = reform(coefv, nx, ny, nzt, /over) 
    122          coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     126         coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    123127         terrev = where(coefv EQ 0) 
    124128         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
    125129 
    126          tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] 
     130         tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 
    127131         div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt) 
    128132         div = reform(div, nx, ny, nzt, /over) 
     
    137141; mise a !values.f_nan de la bordure 
    138142;------------------------------------------------------------ 
    139          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 
    140144            psi(0, *, *) = !values.f_nan 
    141145            psi(nx-1, *, *) = !values.f_nan 
     
    150154; pour le trace graphique 
    151155;------------------------------------------------------------ 
    152          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f'] 
    153          if keyword_set(direc) then psi = moyenne(psi,direc,/nan, boite = boite) 
     156         domdef, (glagmt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
     157         if keyword_set(direc) then psi = moyenne(psi,direc,/nan) 
    154158 
    155159;---------------------------------------------------------------------------- 
     
    170174             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    171175               if nxu NE nx then $ 
    172                 if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]  
     176                if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]  
    173177               IF nxv NE nx THEN $ 
    174                 if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
     178                if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 
    175179               IF nyu NE ny THEN $ 
    176                 if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]  
     180                if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]  
    177181               IF nyv NE ny THEN $ 
    178                 if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
     182                if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 
    179183            END 
    180184            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
     
    191195; calcul du rotationnel 
    192196;---------------------------------------------------------------------------- 
    193          coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt] 
     197         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    194198         terreu = where(coefu EQ 0) 
    195199         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
     
    197201         coefu = reform(coefu, nx, ny, jpt, /over) 
    198202; 
    199          coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt] 
     203         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    200204         terrev = where(coefv EQ 0) 
    201205         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
     
    203207         coefv = reform(coefv, nx, ny, jpt, /over) 
    204208; 
    205          tabf = (fmask())[indice2d+jpi*jpj*premierzt]/(e1f[indice2d]*e2f[indice2d]) 
     209         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 
    206210         tabf = temporary(tabf[*])#replicate(1, jpt) 
    207211         tabf = reform(tabf, nx, ny, jpt, /over) 
     
    217221; mise a !values.f_nan de la bordure 
    218222;------------------------------------------------------------ 
    219          if NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     223         if NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    220224            psi(0, *, *) = !values.f_nan 
    221225            psi(nx-1, *, *) = !values.f_nan 
     
    227231         if terref[0] NE -1 then psi[temporary(terref)] = valmask 
    228232;---------------------------------------------------------------------------- 
    229          domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f'] 
    230          if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan, boite = boite) 
     233         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
     234         if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan) 
    231235;---------------------------------------------------------------------------- 
    232236      END 
     
    255259             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN  
    256260               if nxu NE nx then $ 
    257                 if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]  
     261                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]  
    258262               IF nxv NE nx THEN $ 
    259                 if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
     263                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 
    260264               IF nyu NE ny THEN $ 
    261                 if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]  
     265                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]  
    262266               IF nyv NE ny THEN $ 
    263                 if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
     267                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 
    264268            END 
    265269            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 
     
    273277; calcul du rotationnel 
    274278;------------------------------------------------------------ 
    275          coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt] 
     279         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 
    276280         terreu = where(coefu EQ 0) 
    277281         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 
    278          coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt] 
     282         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 
    279283         terrev = where(coefv EQ 0) 
    280284         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 
    281          tabf = (fmask())[indice2d+jpi*jpj*premierzt]/(e1f[indice2d]*e2f[indice2d]) 
     285         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 
    282286; 
    283287         zu = u*temporary(coefu) 
     
    290294; mise a !values.f_nan de la bordure 
    291295;------------------------------------------------------------ 
    292          if  NOT keyword_set(key_periodique) OR nx NE jpi then begin 
     296         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin 
    293297            psi(0, *) = !values.f_nan 
    294298            psi(nx-1, *) = !values.f_nan 
     
    303307; pour le trace graphique 
    304308;------------------------------------------------------------ 
    305          domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f'] 
    306          if keyword_set(direc) then psi = moyenne(psi,direc,/nan, boite = boite) 
    307  
     309         domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 
     310         if keyword_set(direc) then psi = moyenne(psi,direc,/nan) 
    308311 
    309312      END 
Note: See TracChangeset for help on using the changeset viewer.