Changeset 445


Ignore:
Timestamp:
01/19/11 13:03:43 (13 years ago)
Author:
smasson
Message:

improve the implementation of the partial steps + bugfix in bsf

Location:
trunk/SRC
Files:
5 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/Computation/bsf.pro

    r388 r445  
    8181  un = fitintobox(temporary(un)) 
    8282; 
    83   IF keyword_set(key_partialstep) THEN BEGIN 
    84 ; 2D e3u at the bottom of the ocean 
    85 ; see zgr_zps: e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk)) 
    86     e3u_ps = [[[e3t_ps[firstx:lastx, firsty:lasty]]], [[shift(e3t_ps[firstx:lastx, firsty:lasty], -1, 0)]]] 
    87     e3u_ps = min(temporary(e3u_ps), dimension = 3) 
    88     flagdata = NOT (keyword_set(key_periodic) AND nx EQ jpi) AND total(umsk[nx-1, *, *]) NE 0 
    89 ; level of the bottom of the ocean 
    90     bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 
    91     bottom = 0L > ( long(temporary(bottom)) - 1L ) 
    92 ; the bottom of the ocean in 3D index is: 
    93     bottom = lindgen(nx*ny) + temporary(bottom)*nx*ny 
    94 ; 3D e3t array 
    95     e33 = replicate(1., nx*ny) # e3t[firstz:lastz] 
    96 ; apply e3u_ps to e33 at the bottom of the ocean 
    97     e33[temporary(bottom)] = e3u_ps 
    98 ; 3D e2u array 
    99     e23d = (e2u[firstx:lastx, firsty:lasty])[*] # replicate(1., nzt) 
    100 ; e2u*e3u 
    101     e23 = temporary(e23d) * temporary(e33)  
    102   ENDIF ELSE BEGIN  
    103     e23 = (e2u[firstx:lastx, firsty:lasty])[*] # e3t[firstz:lastz] 
    104   ENDELSE 
     83  e23 = e3u_3d(/e2) 
    10584; 
    10685; mask the array 

  • trunk/SRC/Computation/msf.pro

    r388 r445  
    9595  IF keyword_set(key_partialstep) AND total(msk[*, ny-1, *]) NE 0 THEN flagdata = 1 
    9696; 
    97 ; 
    98   IF keyword_set(key_partialstep) THEN BEGIN 
    99  
    100 ; Rebuild the T-point 3D partial steps array from 2D bottom e3t_ps 
    101 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
    102 ; level of the bottom of the ocean 
    103      bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 
    104      bottom = 0L > ( long(temporary(bottom)) - 1L ) 
    105 ; the bottom of the ocean in 3D index is: 
    106      bottom = lindgen(nx*ny) + temporary(bottom)*nx*ny 
    107 ; 3D e3t array 
    108      e3t_3D = replicate(1, nx*ny) # e3t[firstz:lastz] 
    109 ; apply e3t_ps to e3t_3D at the bottom of the ocean 
    110      e3t_3D[temporary(bottom)] = e3t_ps[firstx:lastx, firsty:lasty] 
    111 ; 3D e3t_3D array 
    112      e3t_3D = reform(e3t_3D, nx, ny, nz, /overwrite) 
    113       
    114 ; Rebuild the V-point 3D partial steps array from T-point 3D e3t_3D array  
    115      tmp = shift(e3t_3D, 0, -1, 0) 
    116      e3v_3D = [ [ (temporary(e3t_3D))[*] ], [ (temporary(tmp))[*] ] ] 
    117      e3v_3D = min(temporary(e3v_3D), dimension = 2) 
    118  
    119 ; 3D e1v array 
    120      e1_3D = (e1v[firstx:lastx, firsty:lasty])[*] # replicate(1., nz) 
    121 ; e1v*e3v 
    122      e13 = temporary(e1_3D) * temporary(e3v_3D)  
    123   ENDIF ELSE BEGIN  
    124      e13 = (e1v[firstx:lastx, firsty:lasty])[*] # e3t[firstz:lastz] 
    125   ENDELSE 
     97  e13 = e3v_3d(/e1)  
    12698; 
    12799; mask the array 

  • trunk/SRC/ToBeReviewed/CALCULS/grossemoyenne.pro

    r442 r445  
    399399    IF keyword_set(key_partialstep) THEN BEGIN 
    400400; the top of the ocean floor is 
    401       IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $ 
    402       ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 
     401      bottom = total(tmask[firstx:lastx, firsty:lasty, *], 3) 
    403402; we suppress columns with only ocean or land 
    404       good = where(bottom NE 0 AND bottom NE nz) 
     403      bottom = long(temporary(bottom)) - firstz - 1L 
     404      good = where(bottom GE 0 AND bottom LT nz) 
    405405; the bottom of the ocean in 3D index is: 
    406       bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny 
     406      bottom = lindgen(nx*ny)+nx*ny*temporary(bottom) 
    407407      IF good[0] NE -1 THEN bottom = bottom[good] $ 
    408408      ELSE bottom = -1 
    … …  
    470470        endif 
    471471      end 
    472       (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 
    473         e33 = replicate(1, 1.*nx*ny)#e3 
    474         e33 = reform(e33, nx, ny, nz, /over) 
    475         IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
    476           IF keyword_set(wdepth) THEN $ 
    477             e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
    478           ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
    479         ENDIF 
     472      (dirx eq 0) and (diry eq 0) and (dirz eq 1) : BEGIN 
     473        e33 = e3_3d() 
    480474        echelle = (temporary(e33)*temporary(mask))[*]#replicate(1, jpt) 
    481475        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
    … …  
    524518      end 
    525519      (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 
    526         e133 = e1[*]#e3 
    527         IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
    528           IF keyword_set(wdepth) THEN $ 
    529             e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
    530           ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
    531         ENDIF 
     520        e133 = e3_3d(/e1) 
    532521        echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt) 
    533522        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
    … …  
    546535      end 
    547536      (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 
    548         e233 = e2[*]#e3 
    549         IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
    550           IF keyword_set(wdepth) THEN $ 
    551             e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
    552           ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
    553         ENDIF 
     537        e233 = e3_3d(/e2) 
    554538        echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt) 
    555539        echelle = reform(echelle, nx, ny, nz, jpt, /over) 
    … …  
    568552      end 
    569553      (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 
    570         e1233 = (e1[*]*e2[*])#e3 
    571         IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
    572           IF keyword_set(wdepth) THEN $ 
    573             e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
    574           ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
    575         ENDIF 
     554        e1233 = e3_3d(/e1, /e2) 
    576555        echelle = (temporary(e1233[*])*temporary(mask[*]))#replicate(1, jpt) 
    577556        echelle = reform(echelle, nx, ny, nz, jpt, /over) 

  • trunk/SRC/ToBeReviewed/CALCULS/moyenne.pro

    r442 r445  
    376376    IF keyword_set(key_partialstep) THEN BEGIN 
    377377; the top of the ocean floor is 
    378       IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $ 
    379       ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 
     378      bottom = total(tmask[firstx:lastx, firsty:lasty, *], 3) 
    380379; we suppress columns with only ocean or land 
    381       good = where(bottom NE 0 AND bottom NE nz) 
     380      bottom = long(temporary(bottom)) - firstz - 1L 
     381      good = where(bottom GE 0 AND bottom LT nz) 
    382382; the bottom of the ocean in 3D index is: 
    383       bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny 
     383      bottom = lindgen(nx*ny)+nx*ny*temporary(bottom) 
    384384      IF good[0] NE -1 THEN bottom = bottom[good] $ 
    385385      ELSE bottom = -1 
    … …  
    447447      end 
    448448      (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 
    449         e33 = replicate(1, 1.*nx*ny)#e3 
    450         e33 = reform(e33, nx, ny, nz, /over) 
     449        e33 = e3_3d() 
    451450        IF keyword_set(sshh) THEN e33[*, *, 0] = e33[*, *, 0] + sshh 
    452         IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
    453           IF keyword_set(wdepth) THEN $ 
    454             e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
    455           ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
    456         ENDIF 
    457451        if keyword_set(integration) then divi = 1 else begin 
    458452          divi = e33*mask 
    … …  
    500494      end 
    501495      (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 
    502         e133 = e1[*]#e3 
    503         e133 = reform(e133, nx, ny, nz, /over) 
     496        e133 = e3_3d(/e1) 
    504497        IF keyword_set(sshh) THEN e133[*, *, 0] = e133[*, *, 0] + sshh*e1 
    505         IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
    506           IF keyword_set(wdepth) THEN $ 
    507             e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
    508           ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
    509         ENDIF 
    510498        if keyword_set(integration) then divi = 1 else BEGIN 
    511499          divi = e133*mask 
    … …  
    525513      end 
    526514      (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 
    527         e233 = e2[*]#e3 
    528         e233 = reform(e233, nx, ny, nz, /over) 
     515        e233 = e3_3d(/e2) 
    529516        IF keyword_set(sshh) THEN e233[*, *, 0] = e233[*, *, 0] + sshh*e2 
    530         IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
    531           IF keyword_set(wdepth) THEN $ 
    532             e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
    533           ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
    534         ENDIF 
    535517        if keyword_set(integration) then divi = 1 else BEGIN 
    536518          divi = e233*mask 
    … …  
    550532      end 
    551533      (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 
    552         e1233 = (e1*e2)[*]#e3 
    553         e1233 = reform(e1233, nx, ny, nz, /over) 
     534        e1233 = e3_3d(/e1, /e2) 
    554535        IF keyword_set(sshh) THEN e1233[*, *, 0] = e1233[*, *, 0] + sshh*e1*e2 
    555         IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 
    556           IF keyword_set(wdepth) THEN $ 
    557             e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 
    558           ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 
    559         ENDIF 
    560536        if keyword_set(integration) then divi = 1 else BEGIN 
    561537          if msknan[0] NE -1 then divi = total(e1233*mask*msknan) $ 
Note: See TracChangeset for help on using the changeset viewer.