Changeset 501 for trunk


Ignore:
Timestamp:
09/20/16 20:15:00 (8 years ago)
Author:
smasson
Message:

set of bugfixes...

Location:
trunk/SRC
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/Grid/computegrid.pro

    r499 r501  
    264264;- 
    265265PRO computegrid, startx, starty, stepxin, stepyin, nxin, nyin $ 
    266                  , XAXIS=xaxis, YAXIS=yaxis, ZAXIS=zaxis, ZZAXIS=zzaxis $ 
    267                  , MASK=mask, GLAMBOUNDARY=glamboundary $ 
    268                  , XMINMESH=xminmesh, XMAXMESH=xmaxmesh $ 
    269                  , YMINMESH=yminmesh, YMAXMESH=ymaxmesh $ 
    270                  , ZMINMESH=zminmesh, ZMAXMESH=zmaxmesh $ 
    271                  , ONEARTH=onearth, PERIODIC=periodic $ 
    272                  , PLAIN=plain, SHIFT=shift, STRIDE=stride $ 
    273                  , YREVERSE=yreverse, ZREVERSE=zreverse  $ 
    274                  , FULLCGRID=fullcgrid, XYINDEX=xyindex $ 
    275                  , UBASE2TBASE=ubase2tbase, VBASE2TBASE=vbase2tbase $ 
    276                  , FBASE2TBASE=fbase2tbase, FFBASE2TBASE=ffbase2tbase $ 
    277                  , STRCALLING=strcalling, ROMSH=romsh, _EXTRA=ex 
     266                 , XAXIS = xaxis, YAXIS = yaxis, ZAXIS = zaxis, ZZAXIS = zzaxis $ 
     267                 , MASK = mask, GLAMBOUNDARY = glamboundary $ 
     268                 , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $ 
     269                 , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $ 
     270                 , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $ 
     271                 , ONEARTH = onearth, PERIODIC = periodic $ 
     272                 , PLAIN = plain, SHIFT = shift, STRIDE = stride $ 
     273                 , YREVERSE = yreverse, ZREVERSE = zreverse  $ 
     274                 , FULLCGRID = fullcgrid, XYINDEX = xyindex $ 
     275                 , UBASE2TBASE = ubase2tbase, VBASE2TBASE = vbase2tbase $ 
     276                 , FBASE2TBASE = fbase2tbase, FFBASE2TBASE = ffbase2tbase $ 
     277                 , STRCALLING = strcalling, ROMSH = romsh, _EXTRA = ex 
    278278; 
    279279  compile_opt idl2, strictarrsubs 
    … …  
    432432    nzz = jpk 
    433433  ENDIF ELSE BEGIN 
    434    jpkglo = long(nz) 
     434    jpkglo = long(nz) 
    435435    IF n_elements(zminmesh) NE 0 THEN izminmesh = long(zminmesh[0]) ELSE izminmesh  = 0l 
    436436    IF n_elements(zmaxmesh) NE 0 THEN izmaxmesh = long(zmaxmesh[0]) ELSE izmaxmesh  = jpkglo-1 
    … …  
    570570; check the periodicity if iyminmesh and iymaxmesh have the default definitions... 
    571571  IF NOT keyword_set(plain) AND key_onearth EQ 1 $ 
    572     AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 AND jpj GE 3 AND jpi GE 2 THEN BEGIN 
     572     AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 AND jpj GE 3 AND jpi GE 2 THEN BEGIN 
    573573 
    574574    CASE 1 OF 
    575575      ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $ 
    576         AND array_equal(gphit[1:*, jpj-1], reverse(gphit[1:*, jpj-3])) EQ 1:BEGIN 
     576         AND array_equal(gphit[1:*, jpj-1], reverse(gphit[1:*, jpj-3])) EQ 1:BEGIN 
    577577; T pivot 
    578578        ymaxmesh = -1 
    … …  
    586586      END 
    587587      ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $ 
    588        AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN 
     588         AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN 
    589589; F pivot 
    590590        ymaxmesh = -1 
    … …  
    604604; check x periodicity... 
    605605;==================================================== 
    606 IF n_elements(periodic) NE 0 THEN forcenoperio = 1 - keyword_set(periodic) 
     606  IF n_elements(periodic) NE 0 THEN forcenoperio = 1 - keyword_set(periodic) 
    607607; check the periodicity if ixminmesh and ixmaxmesh have the default definitions... 
    608608  IF NOT keyword_set(plain) AND NOT keyword_set(forcenoperio) AND key_onearth EQ 1 $ 
    … …  
    610610    CASE 0 OF 
    611611      total((glamt[0, *] - glamt[jpi-2, *]) MOD 360) $ 
    612       + total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN 
     612         + total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN 
    613613        xminmesh = 1 
    614614        xmaxmesh = -1 
    … …  
    685685      IF array_equal(sort(xtest), lindgen(jpi)) EQ 0 THEN BEGIN 
    686686        ras = report(['WARNING: we cannot sort the xaxis with a simple shift...', $ 
    687         'we force key_periodic = 0 and key_shift = 0', $ 
    688         'only horizontal plot may be ok...']) 
     687                      'we force key_periodic = 0 and key_shift = 0', $ 
     688                      'only horizontal plot may be ok...']) 
    689689        key_periodic = 0 
    690690        xnotsorted = 1 
    … …  
    751751      ELSE:BEGIN 
    752752        if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $ 
    753           OR (keyword_set(key_periodic) AND key_irregular) then BEGIN 
     753           OR (keyword_set(key_periodic) AND key_irregular) then BEGIN 
    754754          stepxf = (glamt + 720) MOD 360 
    755755          IF jpj EQ 1 THEN stepxf = reform(stepxf, jpi, jpj, /over) 
    … …  
    758758          stepxf = min(abs(stepxf), dimension = 3) 
    759759          IF NOT keyword_set(key_periodic) THEN $ 
    760             stepxf[jpi-1, *] = stepxf[jpi-2, *] 
     760             stepxf[jpi-1, *] = stepxf[jpi-2, *] 
    761761        ENDIF ELSE BEGIN 
    762762          stepxf = shift(glamt, -1, -1) - glamt 
    763763          IF keyword_set(key_periodic) THEN $ 
    764             stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $ 
     764             stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $ 
    765765          ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *] 
    766766        ENDELSE 
    … …  
    797797        IF jpi GT 1 THEN BEGIN 
    798798          if NOT keyword_set(key_periodic) THEN $ 
    799             stepyf[jpi-1, *] = stepyf[jpi-2, *] 
     799             stepyf[jpi-1, *] = stepyf[jpi-2, *] 
    800800          stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2] 
    801801        ENDIF 
    … …  
    831831      e1t = 0.5*(stepxu+shift(stepxu, 1, 0)) 
    832832      IF NOT keyword_set(key_periodic) THEN $ 
    833         e1t[0, *] = e1t[1, *] 
     833         e1t[0, *] = e1t[1, *] 
    834834    ENDIF ELSE e1t = replicate(stepx, jpi, jpj) 
    835835  ENDIF ELSE e1t = replicate(1b, jpi, jpj) 
    … …  
    857857  IF jpj EQ 1 THEN e2t = reform(e2t, jpi, jpj, /over) 
    858858; 
     859;==================================================== 
     860; def of glamu: defined as the middle of T(i,j) T(i+1,j) 
     861;==================================================== 
     862; 
     863  IF keyword_set(irregular) THEN BEGIN 
     864    glamu = glamt + 0.5 * stepxu 
     865    IF keyword_set(glamboundary) AND key_onearth THEN $ 
     866       glamu = glamboundary[0] > temporary(glamu) < glamboundary[1] 
     867  ENDIF ELSE glamu = glamf 
     868; 
     869  IF jpj EQ 1 THEN glamu = reform(glamu, jpi, jpj, /over) 
     870; 
     871;==================================================== 
     872; def of gphiu: defined as the middle of T(i,j) T(i+1,j) 
     873;==================================================== 
     874; 
     875  IF jpi GT 1 THEN BEGIN 
     876                                ; we must compute stepyu: y distance between T(i+1,j) T(i,j) 
     877    IF keyword_set(key_irregular) THEN BEGIN 
     878      stepyu = shift(gphit, -1, 0) - gphit 
     879      IF NOT keyword_set(key_periodic) THEN $ 
     880         stepyu[jpi-1, *] = stepyu[jpi-2, *] 
     881      gphiu = gphit + 0.5 * stepyu 
     882    ENDIF ELSE gphiu = gphit 
     883  ENDIF ELSE gphiu = gphit 
     884  IF key_onearth THEN gphiu = -90. > gphiu < 90. 
     885; 
     886  IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over) 
     887; 
     888;==================================================== 
     889; def of glamv: defined as the middle of T(i,j) T(i,j+1) 
     890;==================================================== 
     891; 
     892  IF jpj GT 1 THEN BEGIN 
     893                                ; we must compute stepxv: x distance between T(i,j) T(i,j+1) 
     894    IF keyword_set(irregular) THEN BEGIN 
     895      IF keyword_set(key_periodic) THEN BEGIN 
     896        stepxv = (glamt + 720) MOD 360 
     897        stepxv = shift(stepxv, 0, -1) - stepxv 
     898        stepxv = [ [[stepxv]], [[stepxv + 360]], [[stepxv - 360]] ] 
     899        stepxv = min(abs(stepxv), dimension = 3) 
     900      ENDIF ELSE stepxv = shift(glamt, 0, -1) - glamt 
     901      stepxv[*, jpj-1] = stepxv[*, jpj-2] 
     902      glamv = glamt + 0.5 * stepxv 
     903      IF keyword_set(glamboundary) AND key_onearth THEN $ 
     904         glamv = glamboundary[0] > temporary(glamv) < glamboundary[1] 
     905    ENDIF ELSE glamv = glamt 
     906  ENDIF ELSE glamv = glamt 
     907; 
     908;==================================================== 
     909; def of gphiv: defined as the middle of T(i,j) T(i,j+1) 
     910;==================================================== 
     911; 
     912  IF keyword_set(key_irregular) THEN $ 
     913     gphiv = gphit + 0.5 * stepyv $ 
     914  ELSE gphiv = gphif 
     915  IF key_onearth THEN gphiv = -90. > gphiv < 90. 
     916; 
     917  IF jpj EQ 1 THEN gphiv = reform(gphiv, jpi, jpj, /over) 
     918 
    859919;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    860   IF keyword_set(fullcgrid) THEN BEGIN 
    861 ; 
    862 ;==================================================== 
    863 ; def of glamu: defined as the middle of T(i,j) T(i+1,j) 
    864 ;==================================================== 
    865 ; 
    866     IF keyword_set(irregular) THEN BEGIN 
    867       glamu = glamt + 0.5 * stepxu 
    868       IF keyword_set(glamboundary) AND key_onearth THEN $ 
    869         glamu = glamboundary[0] > temporary(glamu) < glamboundary[1] 
    870     ENDIF ELSE glamu = glamf 
    871 ; 
    872     IF jpj EQ 1 THEN glamu = reform(glamu, jpi, jpj, /over) 
    873 ; 
    874 ;==================================================== 
    875 ; def of gphiu: defined as the middle of T(i,j) T(i+1,j) 
    876 ;==================================================== 
    877 ; 
    878     IF jpi GT 1 THEN BEGIN 
    879  ; we must compute stepyu: y distance between T(i+1,j) T(i,j) 
    880       IF keyword_set(key_irregular) THEN BEGIN 
    881        stepyu = shift(gphit, -1, 0) - gphit 
    882         IF NOT keyword_set(key_periodic) THEN $ 
    883           stepyu[jpi-1, *] = stepyu[jpi-2, *] 
    884         gphiu = gphit + 0.5 * stepyu 
    885       ENDIF ELSE gphiu = gphit 
    886     ENDIF ELSE gphiu = gphit 
    887   IF key_onearth THEN gphiu = -90. > gphiu < 90. 
    888 ; 
    889   IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over) 
    890 ; 
    891 ;==================================================== 
    892 ; def of glamv: defined as the middle of T(i,j) T(i,j+1) 
    893 ;==================================================== 
    894 ; 
    895     IF jpj GT 1 THEN BEGIN 
    896  ; we must compute stepxv: x distance between T(i,j) T(i,j+1) 
    897       IF keyword_set(irregular) THEN BEGIN 
    898         IF keyword_set(key_periodic) THEN BEGIN 
    899           stepxv = (glamt + 720) MOD 360 
    900           stepxv = shift(stepxv, 0, -1) - stepxv 
    901           stepxv = [ [[stepxv]], [[stepxv + 360]], [[stepxv - 360]] ] 
    902           stepxv = min(abs(stepxv), dimension = 3) 
    903         ENDIF ELSE stepxv = shift(glamt, 0, -1) - glamt 
    904         stepxv[*, jpj-1] = stepxv[*, jpj-2] 
    905         glamv = glamt + 0.5 * stepxv 
    906         IF keyword_set(glamboundary) AND key_onearth THEN $ 
    907           glamv = glamboundary[0] > temporary(glamv) < glamboundary[1] 
    908       ENDIF ELSE glamv = glamt 
    909     ENDIF ELSE glamv = glamt 
    910 ; 
    911 ;==================================================== 
    912 ; def of gphiv: defined as the middle of T(i,j) T(i,j+1) 
    913 ;==================================================== 
    914 ; 
    915     IF keyword_set(key_irregular) THEN $ 
    916       gphiv = gphit + 0.5 * stepyv $ 
    917     ELSE gphiv = gphif 
    918     IF key_onearth THEN gphiv = -90. > gphiv < 90. 
    919 ; 
    920     IF jpj EQ 1 THEN gphiv = reform(gphiv, jpi, jpj, /over) 
     920  IF keyword_set(fullcgrid) THEN BEGIN     
    921921; 
    922922;==================================================== 
    … …  
    925925; 
    926926    IF jpi GT 1 AND n_elements(stepx) NE 1 THEN $ 
    927       e1u = stepxu ELSE e1u = e1t 
     927       e1u = stepxu ELSE e1u = e1t 
    928928; 
    929929;==================================================== 
    … …  
    10041004; 
    10051005  IF keyword_set(key_onearth)  THEN BEGIN 
    1006     e1t = r * !pi/180. * temporary(e1t) * cos(!pi/180.*gphit) 
    1007     IF keyword_set(fullcgrid) THEN BEGIN 
    1008       e1u = r * !pi/180. * temporary(e1u) * cos(!pi/180.*gphiu) 
    1009       e1v = r * !pi/180. * temporary(e1v) * cos(!pi/180.*gphiv) 
    1010       e1f = r * !pi/180. * temporary(e1f) * cos(!pi/180.*gphif) 
    1011     ENDIF 
     1006    IF keyword_set(key_irregular) THEN BEGIN 
     1007; T 
     1008      e1t = map_npoints(glamu, gphiu, shift(glamu, 1, 0), shift(gphiu, 1, 0), /two_by_two) 
     1009      e1t = reform(e1t, jpi, jpj, /over) 
     1010      IF NOT keyword_set(key_periodic) THEN e1t[0, *] = e1t[1, *] 
     1011      e2t = map_npoints(glamv, gphiv, shift(glamv, 0, 1), shift(gphiv, 0, 1), /two_by_two) 
     1012      e2t = reform(e2t, jpi, jpj, /over) 
     1013      e2t[*, 0] = e2t[*, 1] 
     1014      IF keyword_set(fullcgrid) THEN BEGIN 
     1015; U 
     1016        e1u = map_npoints(glamt, gphit, shift(glamt, -1, 0), shift(gphit, -1, 0), /two_by_two) 
     1017        e1u = reform(e1u, jpi, jpj, /over) 
     1018        IF NOT keyword_set(key_periodic) THEN e1u[jpi-1, *] = e1u[jpi-2, *] 
     1019        e2u = map_npoints(glamf, gphif, shift(glamf, 0, 1), shift(gphif, 0, 1), /two_by_two) 
     1020        e2u = reform(e2u, jpi, jpj, /over) 
     1021        e2u[*, 0] = e2u[*, 1] 
     1022; V 
     1023        e1v = map_npoints(glamf, gphif, shift(glamf, 1, 0), shift(gphif, 1, 0), /two_by_two) 
     1024        e1v = reform(e1v, jpi, jpj, /over) 
     1025        IF NOT keyword_set(key_periodic) THEN e1v[0, *] = e1v[1, *] 
     1026        e2v = map_npoints(glamt, gphit, shift(glamt, 0, -1), shift(gphit, 0, -1), /two_by_two) 
     1027        e2v = reform(e2v, jpi, jpj, /over) 
     1028        e2v[*, jpj-1] = e2v[*, jpj-2] 
     1029; F 
     1030        e1f = map_npoints(glamv, gphiv, shift(glamv, -1, 0), shift(gphiv, -1, 0), /two_by_two) 
     1031        e1f = reform(e1f, jpi, jpj, /over) 
     1032        IF NOT keyword_set(key_periodic) THEN e1f[jpi-1, *] = e1f[jpi-2, *] 
     1033        e2f = map_npoints(glamu, gphiu, shift(glamu, 0, -1), shift(gphiu, 0, -1), /two_by_two) 
     1034        e2f = reform(e2f, jpi, jpj, /over) 
     1035        e2f[*, jpj-1] = e2f[*, jpj-2] 
     1036      ENDIF 
     1037 
     1038    ENDIF ELSE BEGIN  
     1039      e1t = r * !pi/180. * temporary(e1t) * cos(!pi/180.*gphit) 
     1040      IF keyword_set(fullcgrid) THEN BEGIN 
     1041        e1u = r * !pi/180. * temporary(e1u) * cos(!pi/180.*gphiu) 
     1042        e1v = r * !pi/180. * temporary(e1v) * cos(!pi/180.*gphiv) 
     1043        e1f = r * !pi/180. * temporary(e1f) * cos(!pi/180.*gphif) 
     1044      ENDIF 
     1045    ENDELSE 
    10121046  ENDIF 
    10131047; 
    … …  
    10421076; z axis 
    10431077; 
    1044     CASE n_elements(zaxis) OF 
    1045       0:BEGIN 
    1046         gdept = 0. 
    1047         key_zreverse = 0 
    1048       END 
    1049       1:BEGIN 
    1050         gdept = zaxis 
    1051         key_zreverse = 0 
    1052       END 
    1053       ELSE:BEGIN 
    1054         gdept = zaxis[zmin:zmax] 
    1055         IF n_elements(zreverse) EQ 0 THEN BEGIN 
    1056           IF jpk GT 1 THEN BEGIN 
    1057             if gdept[0] GT gdept[1] then key_zreverse = 1 ELSE key_zreverse = 0 
    1058           ENDIF ELSE key_zreverse = 0 
    1059         ENDIF ELSE key_zreverse = zreverse 
    1060         IF keyword_set(key_zreverse) THEN gdept = reverse(gdept) 
    1061       END 
    1062     ENDCASE 
    1063 ; 
    1064     if n_elements(gdept) GT 1 then BEGIN 
    1065       stepz = shift(gdept, -1)-gdept 
    1066       stepz[jpk-1] = stepz[jpk-2] 
    1067       gdepw = 0. > (gdept-stepz/2.) 
    1068     ENDIF ELSE BEGIN 
    1069       stepz = 1. 
    1070       gdepw = gdept 
    1071     ENDELSE 
    1072     IF keyword_set(romsh) THEN gdepw = gdept 
     1078  CASE n_elements(zaxis) OF 
     1079    0:BEGIN 
     1080      gdept = 0. 
     1081      key_zreverse = 0 
     1082    END 
     1083    1:BEGIN 
     1084      gdept = zaxis 
     1085      key_zreverse = 0 
     1086    END 
     1087    ELSE:BEGIN 
     1088      gdept = zaxis[zmin:zmax] 
     1089      IF n_elements(zreverse) EQ 0 THEN BEGIN 
     1090        IF jpk GT 1 THEN BEGIN 
     1091          if gdept[0] GT gdept[1] then key_zreverse = 1 ELSE key_zreverse = 0 
     1092        ENDIF ELSE key_zreverse = 0 
     1093      ENDIF ELSE key_zreverse = zreverse 
     1094      IF keyword_set(key_zreverse) THEN gdept = reverse(gdept) 
     1095    END 
     1096  ENDCASE 
     1097; 
     1098  if n_elements(gdept) GT 1 then BEGIN 
     1099    stepz = shift(gdept, -1)-gdept 
     1100    stepz[jpk-1] = stepz[jpk-2] 
     1101    gdepw = 0. > (gdept-stepz/2.) 
     1102  ENDIF ELSE BEGIN 
     1103    stepz = 1. 
     1104    gdepw = gdept 
     1105  ENDELSE 
     1106  IF keyword_set(romsh) THEN gdepw = gdept 
    10731107; 
    10741108;==================================================== 
    … …  
    10761110;==================================================== 
    10771111; 
    1078     e3t = stepz 
    1079     IF n_elements(stepz) GT 1 THEN BEGIN 
    1080       e3w = 0.5*(stepz+shift(stepz, 1)) 
    1081       e3w[0] = 0.5*e3t[0] 
    1082     ENDIF ELSE e3w = e3t 
     1112  e3t = stepz 
     1113  IF n_elements(stepz) GT 1 THEN BEGIN 
     1114    e3w = 0.5*(stepz+shift(stepz, 1)) 
     1115    e3w[0] = 0.5*e3t[0] 
     1116  ENDIF ELSE e3w = e3t 
    10831117; 
    10841118;==================================================== 
    … …  
    12281262    IF jpi GT 1 THEN BEGIN 
    12291263      if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $ 
    1230         OR (keyword_set(key_periodic) AND key_irregular) then BEGIN 
     1264         OR (keyword_set(key_periodic) AND key_irregular) then BEGIN 
    12311265        stepxf = (glamt + 720) MOD 360 
    12321266        stepxf = shift(stepxf, -1, -1) - stepxf 
    … …  
    12341268        stepxf = min(abs(stepxf), dimension = 3) 
    12351269        IF NOT keyword_set(key_periodic) THEN $ 
    1236           stepxf[jpi-1, *] = stepxf[jpi-2, *] 
     1270           stepxf[jpi-1, *] = stepxf[jpi-2, *] 
    12371271      ENDIF ELSE BEGIN 
    12381272        stepxf = shift(glamt, -1, -1) - glamt 
    12391273        IF keyword_set(key_periodic) THEN $ 
    1240           stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $ 
    1241           ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *] 
     1274           stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $ 
     1275        ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *] 
    12421276      ENDELSE 
    12431277      IF jpj GT 1 THEN BEGIN 
    … …  
    12591293      IF jpi GT 1 THEN BEGIN 
    12601294        if NOT keyword_set(key_periodic) THEN $ 
    1261           stepyf[jpi-1, *] = stepyf[jpi-2, *] 
     1295           stepyf[jpi-1, *] = stepyf[jpi-2, *] 
    12621296        stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2] 
    12631297      ENDIF 
    … …  
    13161350; 
    13171351  IF total(tmask) EQ jpi*jpj*jpk $ 
    1318     AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $ 
     1352     AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $ 
    13191353  ELSE BEGIN 
    13201354; are we using ORCA2 ? 
    … …  
    13571391 
    13581392  ccmeshparameters = {filename:strcalling  $ 
    1359           , glaminfo:float(string(glaminfo, format = '(E11.4)')) $ 
    1360           , gphiinfo:float(string(gphiinfo, format = '(E11.4)')) $ 
    1361           , jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $ 
    1362           , jpi:jpi, jpj:jpj, jpk:jpk $ 
    1363           , ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $ 
    1364           , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $ 
    1365           , izminmesh:izminmesh, izmaxmesh:izmaxmesh $ 
    1366           , key_shift:key_shift, key_periodic:key_periodic $ 
    1367           , key_stride:key_stride, key_gridtype:key_gridtype $ 
    1368           , key_yreverse:key_yreverse, key_zreverse:key_zreverse $ 
    1369           , key_partialstep:key_partialstep, key_onearth:key_onearth} 
     1393                      , glaminfo:float(string(glaminfo, format = '(E11.4)')) $ 
     1394                      , gphiinfo:float(string(gphiinfo, format = '(E11.4)')) $ 
     1395                      , jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $ 
     1396                      , jpi:jpi, jpj:jpj, jpk:jpk $ 
     1397                      , ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $ 
     1398                      , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $ 
     1399                      , izminmesh:izminmesh, izmaxmesh:izmaxmesh $ 
     1400                      , key_shift:key_shift, key_periodic:key_periodic $ 
     1401                      , key_stride:key_stride, key_gridtype:key_gridtype $ 
     1402                      , key_yreverse:key_yreverse, key_zreverse:key_zreverse $ 
     1403                      , key_partialstep:key_partialstep, key_onearth:key_onearth} 
    13701404 
    13711405  ccreadparameters = {funclec_name:'read_ncdf' $ 
    1372           , jpidta:jpidta, jpjdta:jpjdta, jpkdta:jpkdta $ 
    1373           , ixmindta:ixmindta, ixmaxdta:ixmaxdta $ 
    1374           , iymindta:iymindta, iymaxdta:iymaxdta $ 
    1375           , izmindta:izmindta, izmaxdta:izmaxdta} 
     1406                      , jpidta:jpidta, jpjdta:jpjdta, jpkdta:jpkdta $ 
     1407                      , ixmindta:ixmindta, ixmaxdta:ixmaxdta $ 
     1408                      , iymindta:iymindta, iymaxdta:iymaxdta $ 
     1409                      , izmindta:izmindta, izmaxdta:izmaxdta} 
    13761410;------------------------------------------------------------ 
    13771411  IF keyword_set(key_performance) EQ 1 THEN $ 
    1378     print, 'time computegrid', systime(1)-time1 
     1412     print, 'time computegrid', systime(1)-time1 
    13791413;------------------------------------------------------------ 
    13801414  return 

  • trunk/SRC/Grid/ncdf_meshroms.pro

    r493 r501  
    273273    nothing = execute(command) 
    274274  ENDFOR 
    275   d1t = 1.e3*(shift(d1t, -1, 0) - d1t) 
    276   d2t = 1.e3*(shift(d2t, 0, -1) - d2t) 
     275  d1t = shift(d1t, -1, 0) - d1t 
     276  d2t = shift(d2t, 0, -1) - d2t 
    277277  for i = 0, n_elements(namebase2)-1 do begin 
    278278    command = namebase2[i]+'t = '+namebase2[i]+'t[0:jpi-2, 0:jpj-2]' 
    … …  
    306306    nothing = execute(command) 
    307307  ENDFOR 
    308   tmpsave = 2. * 1.e3 * d1u[0, 0:jpj-2] 
    309   d1u = 1.e3*(shift(d1u, -1, 0) - d1u) 
    310   d2u = 1.e3*(shift(d2u, 0, -1) - d2u) 
     308  tmpsave = 2. * d1u[0, 0:jpj-2] 
     309  d1u = shift(d1u, -1, 0) - d1u 
     310  d2u = shift(d2u, 0, -1) - d2u 
    311311  for i = 0, n_elements(namebase2)-1 do begin 
    312312    command = namebase2[i]+'u = '+namebase2[i]+'u[*, 0:jpj-2]' 
    … …  
    342342    nothing = execute(command) 
    343343  ENDFOR 
    344   d1v = 1.e3*(shift(d1v, -1, 0) - d1v) 
    345   tmpsave = 2. * 1.e3 * d2v[0:jpi-2, 0] 
    346   d2v = 1.e3*(shift(d2v, 0, -1) - d2v) 
     344  d1v = shift(d1v, -1, 0) - d1v 
     345  tmpsave = 2. * d2v[0:jpi-2, 0] 
     346  d2v = shift(d2v, 0, -1) - d2v 
    347347  for i = 0, n_elements(namebase2)-1 do begin 
    348348    command = namebase2[i]+'v = '+namebase2[i]+'v[0:jpi-2, *]' 
    … …  
    372372    nothing = execute(command) 
    373373  ENDFOR 
    374   tmpsave1 = 2. * 1.e3 * d1f[0, *] 
    375   d1f = 1.e3*(shift(d1f, -1, 0) - d1f) 
    376   tmpsave2 = 2. * 1.e3 * d2f[*, 0] 
    377   d2f = 1.e3*(shift(d2f, 0, -1) - d2f) 
     374  tmpsave1 = 2. * d1f[0, *] 
     375  d1f = shift(d1f, -1, 0) - d1f 
     376  tmpsave2 = 2. * d2f[*, 0] 
     377  d2f = shift(d2f, 0, -1) - d2f 
    378378  fmaskredy = byte(maskf[jpi-1, *]) 
    379379  IF jpk GT 1 THEN fmaskredy = reform(fmaskredy[*]#replicate(1b, jpk), 1, jpj, jpk, /overwrite) 

  • trunk/SRC/Grid/romsdepth.pro

    r370 r501  
    2525; 
    2626;- 
    27 FUNCTION romsdepth 
     27FUNCTION romsdepth, new_s_coord = new_s_coord 
    2828; 
    2929  compile_opt idl2, strictarrsubs 
    … …  
    4646  grille, -1, -1, -1, -1, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz 
    4747  hroms = hroms[firstx:lastx, firsty:lasty] 
    48   nt = n_elements(zeta)/nx/ny 
    49 ; 
    50   cff1 = 1./sinh(theta_s) 
    51   cff2 = 0.5/tanh(0.5*theta_s) 
     48  nt = n_elements(zeta[firstx:lastx, firsty:lasty, *])/nx/ny 
    5249; 
    5350  IF type EQ 'W' THEN BEGIN 
    … …  
    5956  ENDELSE 
    6057; 
    61   cs = (1.-theta_b)*cff1*sinh(theta_s*sc)+theta_b*(cff2*tanh(theta_s*(sc+0.5))-0.5) 
    62   cff = hc*(sc-cs) 
    63   cff1 = cs 
     58  IF keyword_set(new_s_coord) THEN BEGIN 
     59 
     60    IF theta_s GT 0. THEN BEGIN 
     61      csrf = ( 1. - cosh(theta_s*sc) ) / ( cosh(theta_s) - 1. ) 
     62    ENDIF ELSE BEGIN 
     63      csrf = - sc*sc 
     64    ENDELSE  
     65    if theta_b GT 0. THEN BEGIN  
     66      cs = ( exp(theta_b*csrf) - 1. ) / ( 1. - exp(-theta_b) )   
     67    ENDIF ELSE BEGIN 
     68      cs = csrf 
     69    ENDELSE 
     70     
     71    hh = (hroms[*]#replicate(1., jpk)) 
     72    z0 = hh * ( replicate(1., nx*ny)#(hc*sc)[*] + hh * replicate(1., nx*ny)#cs[*] ) / ( hc + hh ) 
     73     
     74     
     75  ENDIF ELSE BEGIN 
     76     
     77    cff1 = 1./sinh(theta_s) 
     78    cff2 = 0.5/tanh(0.5*theta_s) 
     79 
     80    cs = (1.-theta_b)*cff1*sinh(theta_s*sc)+theta_b*(cff2*tanh(theta_s*(sc+0.5))-0.5) 
     81    cff = hc*(sc-cs) 
     82    cff1 = cs 
     83 
     84    z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk)) 
     85 
     86  ENDELSE 
     87 
    6488; 
    6589  hinv = 1./hroms 
    6690  hinv = hinv[*]#replicate(1., jpk) 
    6791; put a z dimension to zeta 
    68   zeta = transpose(temporary(zeta)) 
     92  zeta = transpose((temporary(zeta))[firstx:lastx, firsty:lasty, *]) 
    6993  zeta = reform((temporary(zeta))[*]#replicate(1., jpk), nt, ny, nx, jpk, /overwrite) 
    7094  zeta = transpose(temporary(zeta), [2, 1, 3, 0]) 
    7195 
    72   z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk)) 
    7396  z = (z0[*]#replicate(1., nt)) + temporary(zeta) * ((1.+z0*hinv)[*]#replicate(1., nt)) 
    7497  z = reform(z, nx, ny, jpk, nt, /overwrite) 

  • trunk/SRC/Interpolation/compute_fromreg_imoms3_weigaddr.pro

    r495 r501  
    88; Interpolation 
    99; 
    10 ; @param alonin {in}{required}{type=2d array} 
     10; @param alonin {in}{required}{type=1d array} 
    1111; longitude of the input data 
    1212; 
    13 ; @param alatin {in}{required}{type=2d array} 
     13; @param alatin {in}{required}{type=1d array} 
    1414; latitude of the input data 
    1515; 
    … …  
    133133  IF total(alon[indexlon] GT olon) NE 0 THEN stop 
    134134  IF total(alon[indexlon + 1L] LE olon) NE 0 THEN stop 
    135   IF (where(indexlon LE 1L     ))[0] NE -1 THEN stop 
     135  IF (where(indexlon LT 1L     ))[0] NE -1 THEN stop 
    136136  IF (where(indexlon GE jpia-3L))[0] NE -1 THEN stop 
    137137  indexlat = value_locate(alat, olat) 

  • trunk/SRC/Interpolation/neighbor.pro

    r371 r501  
    6969    IF arg_present(distance) THEN distance = sqrt(distance) 
    7070  ENDELSE 
    71   RETURN, where(distance EQ min(distance)) 
     71  minval = min(distance, indmin) 
     72  RETURN, indmin 
    7273END 

  • trunk/SRC/ReadWrite/ncdf_getatt.pro

    r495 r501  
    5858; if this attribute does not exist. 
    5959; 
     60; @keyword interval_operation 
     61; Set this keyword to a named variable in which the value of 
     62; interval_operation is returned. Return -1 if this attribute 
     63; does not exist. 
     64; 
    6065; @keyword _EXTRA 
    6166; defined only to be able to call <pro>ncdf_getatt</pro> with the 
    … …  
    7681                 , MISSING_VALUE = missing_value, UNITS = units, CALENDAR = calendar $ 
    7782                 , LONG_NAME = long_name, STANDARD_NAME = standard_name, DOUBLE = double $ 
    78                  , _extra = ex 
     83                 , INTERVAL_OPERATION = interval_operation, _extra = ex 
    7984; 
    8085  compile_opt idl2, strictarrsubs 
    … …  
    9499  long_name = '' 
    95100  standard_name = '' 
     101  interval_operation = -1 
    96102; 
    97103  varinq = ncdf_varinq(cdfid, varid) 
    … …  
    136142        '_fillvalue':ncdf_attget, cdfid, varid, attname, missing_value 
    137143        'missing_value':ncdf_attget, cdfid, varid, attname, missing_value 
     144        'interval_operation':BEGIN 
     145          ncdf_attget, cdfid, varid, attname, tmp 
     146          tmp = string(tmp) 
     147          interval_operation = fix(strmid(tmp,0,strpos(tmp,' '))) 
     148        END 
    138149        ELSE: 
    139150      ENDCASE 
    … …  
    148159      missing_value = tmp 
    149160    ENDIF ELSE BEGIN 
    150       ras = report('Warning: missing value is not a number: ' + missing_value) 
    151       missing_value = 'no' 
     161      ;; ras = report('Warning: missing value is not a number: ' + missing_value) 
     162      ;; missing_value = 'no' 
     163      missing_value = (byte(strlowcase(a)))[0] 
    152164    ENDELSE 
    153165  ENDIF 

  • trunk/SRC/ToBeReviewed/GRILLE/grille.pro

    r495 r501  
    396396; for the vertical sections with roms 
    397397  IF keyword_set(type) AND n_elements(romszinfos) NE 0 THEN BEGIN 
    398     romsdp = romsdepth() 
     398    romsdp = romsdepth(/new_s_coord) 
    399399    IF romsdp[0] NE -1 THEN BEGIN 
    400400      IF jpt EQ 1 THEN BEGIN 

  • trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro

    r500 r501  
    120120                    , GRID = grid, CALLITSELF = callitself, DIREC = direc $ 
    121121                    , ZETAFILENAME = zetafilename, ZETAZERO = zetazero $ 
    122                     , ZINVAR = zinvar, _EXTRA = ex 
     122                    , ZINVAR = zinvar, ISROMS = isroms, _EXTRA = ex 
    123123; 
    124124  compile_opt idl2, strictarrsubs 
    … …  
    164164    dimnames[i] = tmp 
    165165  ENDFOR 
    166 ;------------------------------------------------------------ 
    167 ; we check if the variable has a vertical dimension and/or 
    168 ; a record dimension. This is useful to define boxzoom 
    169 ; keyword in domdef 
    170 ;------------------------------------------------------------ 
    171   dummy = where(varinq.dim EQ inq.recdim, tinvar) 
    172 ; check the presence of a vertical dimension according to the 
    173 ; number of dimensions and the presence of a record dimension 
    174   zinvar = (varinq.ndims EQ 3 AND tinvar NE 1) OR varinq.ndims EQ 4 
    175 ;------------------------------------------------------------ 
    176 ; shall we redefine the grid parameters 
    177 ;------------------------------------------------------------ 
    178   IF keyword_set(init) THEN initncdf, filename, _extra = ex 
    179166;------------------------------------------------------------ 
    180167; check the time axis and the debut and ending dates 
    … …  
    248235    END 
    249236  ENDCASE 
     237  IF inq.recdim EQ -1 AND jpt NE 1 THEN inq.recdim = varinq.dim[varinq.ndims-1] 
     238;------------------------------------------------------------ 
     239; we check if the variable has a vertical dimension and/or 
     240; a record dimension. This is useful to define boxzoom 
     241; keyword in domdef 
     242;------------------------------------------------------------ 
     243  dummy = where(varinq.dim EQ inq.recdim, tinvar) 
     244; check the presence of a vertical dimension according to the 
     245; number of dimensions and the presence of a record dimension 
     246  zinvar = (varinq.ndims EQ 3 AND tinvar NE 1) OR varinq.ndims EQ 4 
     247;------------------------------------------------------------ 
     248; shall we redefine the grid parameters 
     249;------------------------------------------------------------ 
     250  IF keyword_set(init) THEN initncdf, filename, _extra = ex 
    250251;------------------------------------------------------------ 
    251252; Name of the grid on which the field refer to. 
    … …  
    405406; if it is roms outputs, we need to get additional infos... 
    406407  IF NOT keyword_set(callitself) THEN BEGIN 
    407     IF strmid(dimnames[0], 0, 3) EQ 'xi_' AND strmid(dimnames[1], 0, 4) EQ 'eta_' THEN BEGIN 
     408    IF (strmid(dimnames[0], 0, 3) EQ 'xi_' AND strmid(dimnames[1], 0, 4) EQ 'eta_') OR keyword_set(isroms) THEN BEGIN 
    408409      ncdf_attget, cdfid, 'theta_s', theta_s, /global 
    409410      ncdf_attget, cdfid, 'theta_b', theta_b, /global 

  • trunk/SRC/ToBeReviewed/PLOTS/DESSINE/sbar_plot.pro

    r495 r501  
    5353   if NOT keyword_set(NOREINITPLT) then reinitplt, _extra = ex 
    5454; 2) I place the drawing on the screen like on the postscript 
    55    IF chkstru(ex, 'overplot') EQ 0 THEN placedessin, 'autre', _extra = ex $ 
    56    ELSE IF ex.overplot EQ 0 THEN placedessin, 'autre', _extra = ex 
     55   IF chkstru(ex, 'overplot') EQ 0 THEN placedessin, 'autre', /nocolorbar, _extra = ex $ 
     56   ELSE IF ex.overplot EQ 0 THEN placedessin, 'autre', /nocolorbar, _extra = ex 
    5757; 3) Drawing 
    5858   if n_elements(COLORS) NE 0 then BEGIN 

  • trunk/SRC/ToBeReviewed/PLOTS/DESSINE/splot.pro

    r495 r501  
    4545   if NOT keyword_set(NOREINITPLT) then reinitplt, _extra = ex 
    4646; 2) i put the drawing on the screen like on the postscript 
    47    placedessin, 'autre', _extra = ex 
     47   placedessin, 'autre', /nocolorbar, _extra = ex 
    4848; 3) Drawing 
    4949   if n_elements(y) EQ 0 then plot,  x, xstyle = 1, ystyle = 1, _EXTRA = ex $ 
Note: See TracChangeset for help on using the changeset viewer.