Gaelle 2/12/2015 : MAJ 004_Reunion
authorGaelle Tanguy <gaelle.tanguy@meteo.fr>
Wed, 2 Dec 2015 15:32:26 +0000 (15:32 +0000)
committerPhilippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Thu, 19 May 2016 14:44:49 +0000 (16:44 +0200)
24 files changed:
MY_RUN/KTEST/004_Reunion/001_prep_pgd/PRE_PGD1.nam
MY_RUN/KTEST/004_Reunion/001_prep_pgd/clean_prep_pgd_xyz
MY_RUN/KTEST/004_Reunion/002_prep_ideal_case/PRE_IDEA1.nam
MY_RUN/KTEST/004_Reunion/002_prep_ideal_case/clean_prep_ideal_case_xyz
MY_RUN/KTEST/004_Reunion/002_prep_ideal_case/run_prep_ideal_case_xyz
MY_RUN/KTEST/004_Reunion/003_mesonh/EXSEG1.nam
MY_RUN/KTEST/004_Reunion/003_mesonh/EXSEG1.nam_CEN4TH
MY_RUN/KTEST/004_Reunion/003_mesonh/EXSEG1.nam_WENO
MY_RUN/KTEST/004_Reunion/003_mesonh/run_mesonh_xyz
MY_RUN/KTEST/004_Reunion/005_diaprog/dir_series [deleted file]
MY_RUN/KTEST/004_Reunion/005_diaprog/run_diaprog
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/MESONHtools.ncl [new file with mode: 0644]
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/clean_ncl [new file with mode: 0755]
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/get_ncl_highres_files [new file with mode: 0755]
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_BasicMap.ncl [new file with mode: 0644]
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_Cloud.ncl [new file with mode: 0644]
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_CrossSection.ncl [new file with mode: 0644]
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_ModelLevels.ncl [new file with mode: 0644]
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_PressureLevel.ncl [new file with mode: 0644]
MY_RUN/KTEST/004_Reunion/005_ncl_nc4/run_ncl [new file with mode: 0755]
MY_RUN/KTEST/004_Reunion/006_ncl/clean_ncl [new file with mode: 0755]
MY_RUN/KTEST/004_Reunion/006_ncl/plot_Reunion.ncl [new file with mode: 0644]
MY_RUN/KTEST/004_Reunion/006_ncl/run_ncl [new file with mode: 0755]
MY_RUN/KTEST/004_Reunion/Makefile

index 1316b5c..eace45f 100644 (file)
@@ -5,7 +5,7 @@
 &NAM_CONFZ 
  ! NZ_VERB=5 , NB_PROCIO_W=8 
 /
-&NAM_PGDFILE CPGDFILE='REUNION_PGD_1km5_520' /
+&NAM_PGDFILE CPGDFILE='REUNION_PGD_1km5' /
 &NAM_CONF_PROJ
  XLAT0=-20., XLON0=55., XRPK=0., XBETA=0. /
 &NAM_CONF_PROJ_GRID
index c42b280..69644ac 100755 (executable)
@@ -1,2 +1,2 @@
 set -x
-rm -f  *.dir *.hdr REUNION_PGD_1km5* OUTPUT_* pipe* *.tex
+rm -f  *.dir *.hdr REUNION_PGD_1km5* OUTPUT_* pipe* *.tex file_for_xtransfer
index fecdb79..a30e8bb 100644 (file)
@@ -2,7 +2,7 @@
 &NAM_CONFZ 
   ! NZ_VERB=5 , NB_PROCIO_R=8 , NB_PROCIO_W=8 
 /
-&NAM_REAL_PGD CPGD_FILE = 'REUNION_PGD_1km5_520',
+&NAM_REAL_PGD CPGD_FILE = 'REUNION_PGD_1km5',
               LREAD_ZS= T, LREAD_GROUND_PARAM= T /
 &NAM_CONF_PRE LCARTESIAN= F
               LBOUSS= F  CEQNSYS='LHE'
@@ -11,7 +11,7 @@
               !JPHEXT = 3 , NHALO = 3 
               /
 &NAM_CONFn LUSERV=T /
-&NAM_LUNITn CINIFILE = 'REUNION_IDEA_520',CINIFILEPGD = 'REUNION_PGD_1km5_520' /
+&NAM_LUNITn CINIFILE = 'REUNION_IDEA',CINIFILEPGD = 'REUNION_PGD_1km5' /
 &NAM_DYNn_PRE 
   ! CPRESOPT='ZRESI' 
 /
index 60972da..a87860c 100755 (executable)
@@ -1,2 +1,2 @@
 set -x
-rm -f REUNION_PGD_1km5_410* REUNION_IDEA_410* OUTPUT_LISTING* pipe* *.tex
+rm -f REUNION_PGD_1km5* REUNION_IDEA* OUTPUT_LISTING* pipe* *.tex file_for_xtransfer
index d8c75bd..be2dd65 100755 (executable)
@@ -4,8 +4,8 @@
 #MNH_LIC for details. version 1.
 set -x
 set -e
-ln -sf ../001_prep_pgd/REUNION_PGD_1km5_520.*??? .
-rm -f REUNION_IDEA_520.* OUTPUT_LISTING* pipe* *.tex
+ln -sf ../001_prep_pgd/REUNION_PGD_1km5.*??? .
+rm -f REUNION_IDEA* OUTPUT_LISTING* pipe* *.tex
 
 time ${MONORUN} PREP_IDEAL_CASE${XYZ}          
 #ddd --directory=~/../src/dir_obj PREP_IDEAL_CASE${XYZ}
index 91cd055..b3d90dd 100644 (file)
@@ -1,15 +1,15 @@
 &NAM_CONFIO  LCDF4=T, LLFIOUT=T, LLFIREAD=F /
 &NAM_CONFZ
 ! NZ_VERB=5 , NB_PROCIO_R=8 , NB_PROCIO_W=8
+ ! NZ_VERB=5 , NB_PROCIO_R=8 , NB_PROCIO_W=8
 /
-&NAM_LUNITn  CINIFILE = "REUNION_IDEA_520",
-             CINIFILEPGD="REUNION_PGD_1km5_520" /
+&NAM_LUNITn  CINIFILE = "REUNION_IDEA",
+             CINIFILEPGD="REUNION_PGD_1km5" /
 &NAM_DYNn  XTSTEP = 20.,  
            !CPRESOPT = "ZRESI", 
            LITRADJ = T,
            LHORELAX_UVWTH = T, LHORELAX_RV = T, LVE_RELAX = T,
            NRIMX = 5, NRIMY = 5, XRIMKMAX = 0.01, XT4DIFU = 500. /
-&NAM_ADVn  CUVW_ADV_SCHEME = "WENO_K" , CMET_ADV_SCHEME = "PPM_01" /
+&NAM_ADVn  CUVW_ADV_SCHEME = "WENO_K" , CMET_ADV_SCHEME = "PPM_01",NWENO_ORDER=5,CTEMP_SCHEME='RK53' /
 &NAM_PARAMn  CTURB = "TKEL", CRAD = "FIXE",
              CCLOUD = "NONE", CDCONV = "NONE" /
 &NAM_LBCn  CLBCX = 2*"OPEN", CLBCY = 2*"OPEN" /
index cf5d50f..0f220a4 100644 (file)
@@ -1,8 +1,9 @@
+&NAM_CONFIO  LCDF4=T, LLFIOUT=T, LLFIREAD=F /
 &NAM_CONFZ
- ! NZ_VERB=5 , NB_PROCIO_R=8 , NB_PROCIO_W=1
+ ! NZ_VERB=5 , NB_PROCIO_R=8 , NB_PROCIO_W=8
 /
-&NAM_LUNITn  CINIFILE = "REUNION_IDEA_410",
-             CINIFILEPGD="REUNION_PGD_1km5_410" /
+&NAM_LUNITn  CINIFILE = "REUNION_IDEA",
+             CINIFILEPGD="REUNION_PGD_1km5" /
 &NAM_DYNn  XTSTEP = 5.,  
            !CPRESOPT = "ZRESI", 
            LITRADJ = T,
@@ -18,7 +19,9 @@
              NBJSLICE=1 NJSLICEL(1)=30 NJSLICEH(1)=35 /
 &NAM_CONF  CCONF = "START", NMODEL = 1,
            NVERB = 5, CEXP = "REUNI", CSEG = "00A20" ,
-           CSPLIT="BSPLITTING" /
+           CSPLIT="BSPLITTING" 
+           !JPHEXT =3 , NHALO = 3 
+           /
 &NAM_DYN  XSEGLEN = 40.,
           LCORIO = F, LNUMDIFU = T,
           XALKTOP = 0.01, XALZBOT = 14000. /
index f588e40..d3f6d77 100644 (file)
@@ -1,14 +1,15 @@
+&NAM_CONFIO  LCDF4=T, LLFIOUT=T, LLFIREAD=F /
 &NAM_CONFZ
- ! NZ_VERB=5 , NB_PROCIO_R=8 , NB_PROCIO_W=1
+ ! NZ_VERB=5 , NB_PROCIO_R=8 , NB_PROCIO_W=8
 /
-&NAM_LUNITn  CINIFILE = "REUNION_IDEA_410",
-             CINIFILEPGD="REUNION_PGD_1km5_410" /
+&NAM_LUNITn  CINIFILE = "REUNION_IDEA",
+             CINIFILEPGD="REUNION_PGD_1km5" /
 &NAM_DYNn  XTSTEP = 20.,  
            !CPRESOPT = "ZRESI", 
            LITRADJ = T,
            LHORELAX_UVWTH = T, LHORELAX_RV = T, LVE_RELAX = T,
            NRIMX = 5, NRIMY = 5, XRIMKMAX = 0.01, XT4DIFU = 500. /
-&NAM_ADVn  CUVW_ADV_SCHEME = "WENO_K" , CMET_ADV_SCHEME = "PPM_01" /
+&NAM_ADVn  CUVW_ADV_SCHEME = "WENO_K" , CMET_ADV_SCHEME = "PPM_01", NWENO_ORDER=5,CTEMP_SCHEME='RK53'/
 &NAM_PARAMn  CTURB = "TKEL", CRAD = "FIXE",
              CCLOUD = "NONE", CDCONV = "NONE" /
 &NAM_LBCn  CLBCX = 2*"OPEN", CLBCY = 2*"OPEN" /
@@ -18,7 +19,9 @@
              NBJSLICE=1 NJSLICEL(1)=30 NJSLICEH(1)=35 /
 &NAM_CONF  CCONF = "START", NMODEL = 1,
            NVERB = 5, CEXP = "REUNI", CSEG = "00A20" ,
-           CSPLIT="BSPLITTING" /
+           CSPLIT="BSPLITTING" 
+           !JPHEXT =3 , NHALO = 3 
+           /
 &NAM_DYN  XSEGLEN = 40.,
           LCORIO = F, LNUMDIFU = F,
           XALKTOP = 0.01, XALZBOT = 14000. /
index 910473e..1b94aa9 100755 (executable)
@@ -4,7 +4,7 @@
 #MNH_LIC for details. version 1.
 set -x
 set -e
-ln -fs ../002_prep_ideal_case/REUNION_IDEA_520.*??? .
-ln -sf ../001_prep_pgd/REUNION_PGD_1km5_520.* .
+ln -fs ../002_prep_ideal_case/REUNION_IDEA.*??? .
+ln -sf ../001_prep_pgd/REUNION_PGD_1km5.* .
 rm -f REUNI.1.* OUT*
 time ${MPIRUN} MESONH${XYZ}
diff --git a/MY_RUN/KTEST/004_Reunion/005_diaprog/dir_series b/MY_RUN/KTEST/004_Reunion/005_diaprog/dir_series
deleted file mode 100644 (file)
index f3ac63a..0000000
+++ /dev/null
@@ -1,63 +0,0 @@
-LINVWB=T
-LCOLAREA=T LISO=F
-XLWCONT=2
-LMINMAX=T
-! label de toutes les iso
-NULBLL=0
-! pour obtenir des labels d isolignes
-XSIZEL=0.012
-! localisation du minmax dans le print minmax
-LMNMXLOC=T
-! High&Low
-NHI=0
-! label des axes en indices plutot que km en CH
-LINDAX=T
-NIGRNC=10
-! vecteur
-XVRL=0.15
-NISKIP=4
-! pas de dilatation de la composante W
-LDILW=F
-! min max vecteur
-LVECTMNMX=T
-! min max du zoom si NIMNMX<=0
-LCVZOOM=T
-! normaliser la longueur des vecteurs a 20m/s et tracer tous les vecteurs
-XVHC=-20
-LVSUPSCA=T
-! pour agrandir la legende a droite
-LVPTUSER=T
-XVPTL=0.1 XVPTR=0.9 !XVPTT=0.8
-!LVPTVUSER=T
-!XVPTVL=0.08 XVPTVR=0.78 XVPTVT=0.8
-!LVPTPVUSER=T
-!XVPTPVL=0.08 XVPTPVR=0.78 XVPTPVT=0.8
-!LVPTXYUSER=T
-!XVPTXYL=0.08 XVPTXYR=0.78 XVPTXYT=0.8
-   
-_file_'REUNI.1.00A20.000'
-visu
-print group
-!
-print TSERIES dim proc time
-!print TSERIES minmax
-TSERIES_FT1_
-!print TSERIES(1:1,1:1,1:1,1:35,1:1,1:1)
-!print TSERIES(1:1,1:1,1:1,1:35,1:1,2:2)
-!print TSERIES(1:1,1:1,1:1,1:35,1:1,3:3)
-!print TSERIES(1:1,1:1,1:1,1:35,1:1,9:9)
-!
-print ZTSERIES dim proc
-!ZTSERIES_PVT_
-!
-print XTSERIES01 dim proc
-print XTSERIES02 dim proc
-print XTSERIES03 dim proc
-!XTSERIES01_PXT_
-!XTSERIES02_PXT_
-!XTSERIES03_PXT_
-!XTSERIES04_PXT_
-!XTSERIES05_PXT_
-! .... ?!
-!XTSERIES05_P_PROC1_PYT_
-quit
index ec9b9cd..e925a17 100755 (executable)
@@ -8,6 +8,3 @@ ln -sf ../004_convdia/REUNI.1.00A20.004dg.??? .
 rm -f dir.0*
 ${POSTRUN} diaprog < dir_Reunion
 mv gmeta gmeta_Reunion || echo pas de gmeta
-ln -sf ../003_mesonh/REUNI.1.00A20.000.??? .
-${POSTRUN} diaprog < dir_series
-mv gmeta gmeta_series || echo pas de gmeta
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/MESONHtools.ncl b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/MESONHtools.ncl
new file mode 100644 (file)
index 0000000..6f810bb
--- /dev/null
@@ -0,0 +1,915 @@
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"\r
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"\r
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"\r
+load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"\r
+\r
+;-------------------------------------------------------------\r
+;contains:\r
+; procedure MESONH_map_c\r
+;function mnh_map(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical)   \r
+;function mnh_map_overlays(in_file[1]:file,wks:graphic,plots[*]:graphic, \\r
+;                                  opt_arg[1]:logical,opt_mp[1]:logical)\r
+;function MESONH_pinter( pfield:numeric, loc_param:numeric, ppabs:numeric )\r
+;-------------------------------------------------------------\r
+\r
+;==============================================================\r
+; J.-P. CHABOUREAU \r
+; This is a driver that selects the appropriate \r
+; mapping function based upon the file variables RPK, BETA, LATOR, LONOR\r
+;\r
+;\r
+; Sample usage:\r
+;             a = addfile("...", r")\r
+;             IMAX = a->IMAX\r
+;             JMAX = a->JMAX\r
+;             lat2d = new((/JMAX,IMAX/),"double")\r
+;             lat2d(:,:)=0.\r
+;             lon2d = new((/JMAX,IMAX/),"double")\r
+;             lon2d(:,:)=0.\r
+;             icorners = new((/2,2/),"integer")\r
+;             icorners(:,:)=0\r
+;             res  = True\r
+;             MESONH_map_c (a, res, lat2d, lon2d, icorners)\r
+;\r
+;\r
+undef("MESONH_map_c")\r
+;==============================================================\r
+procedure MESONH_map_c (in_file:file, res:logical, plat, plon, icorner)\r
+;==============================================================\r
+;local rank, dimll, nlat, mlon, lat, lon\r
+local rank, dimll, nlat, mlon\r
+begin\r
+\r
+; Check if the variable RPK is in the file\r
+; ----------------------------------------\r
+if(isfilevar(in_file,"RPK")) then\r
+\r
+; Read projection parameters\r
+; -------------------------\r
+   ZRPK   = in_file->RPK\r
+   ZLATOR = in_file->LATOR\r
+   ZLONOR = in_file->LONOR\r
+   ZBETA = in_file->BETA\r
+   ZLAT0 = in_file->LAT0\r
+   ZLON0 = in_file->LON0\r
+\r
+; Case netcdf from lfi2cdf\r
+; -------------------------\r
+\r
+   if(isfilevar(in_file,"IMAX"))\r
+       XHAT=in_file->XHAT\r
+       YHAT=in_file->YHAT\r
+       IMAX= dimsizes(XHAT)-2\r
+       JMAX= dimsizes(YHAT)-2\r
+       zdx=XHAT(2)-XHAT(1)\r
+       zdy=YHAT(2)-YHAT(1)\r
+\r
+; unstagger\r
+       do ji=0,IMAX-1\r
+         XHAT(ji)=XHAT(ji)+zdx*1.5\r
+       end do\r
+       do jj=0,JMAX-1\r
+         YHAT(jj)=YHAT(jj)+zdy*1.5\r
+       end do\r
+\r
+   else\r
+\r
+; Case netcdf from extractdia\r
+; ---------------------------\r
+       XHAT=in_file->W_E_direction\r
+       YHAT=in_file->S_N_direction\r
+       IMAX= dimsizes(XHAT)\r
+       JMAX= dimsizes(YHAT)\r
+       zdx=XHAT(2)-XHAT(1)\r
+       zdy=YHAT(2)-YHAT(1)\r
+\r
+   end if\r
+\r
+  print ("LATOR="+ZLATOR+" - LONOR="+ZLONOR)\r
+  print ("ZLAT0="+ZLAT0+" - ZLON0="+ZLON0)\r
+  print ("ZDX="+zdx+" - RPK="+ZRPK+" - BETA="+ZBETA)\r
+  print ("IMAX="+IMAX+" - JMAX="+JMAX)\r
+\r
+  if (ZRPK.gt.0) \r
+  ;   Stereographic projection\r
+; ---------------------------\r
+    res@mpProjection = "Stereographic"\r
+    res@mpCenterLonF           = ZLON0\r
+    res@mpCenterRotF = ZBETA\r
+    res@mpCenterLatF           = 90. \r
+  end if\r
+\r
+  if (ZRPK.lt.0) \r
+  ;   Stereographic projection\r
+; ---------------------------\r
+    res@mpProjection = "Stereographic"\r
+    res@mpCenterLonF           = ZLON0\r
+    res@mpCenterRotF = ZBETA\r
+    res@mpCenterLatF           = -90. \r
+  end if\r
+\r
+  if (ZRPK.eq.0) then\r
+  ;   Mercator projection\r
+; ---------------------------\r
+    res@mpProjection = "Mercator"\r
+  end if\r
+\r
+ print("Map projection="+res@mpProjection)\r
+\r
+else \r
+  print ("MESONH_map_c: Error no RPK variable in input file")\r
+end if\r
+\r
+;=================================================;\r
+; calculate 2D lat and lon \r
+; based on src/mesonh_MOD/mode_gridproj.f90\r
+;=================================================;\r
+\r
+; Constants  \r
+; -----------\r
+  if(isfilevar(in_file,"IMAX"))\r
+  XRADIUS=6371229.0d ; Earth radius (meters)\r
+  else\r
+  XRADIUS=6371.2290d ; Earth radius (km)\r
+  end if\r
+  XPI=2.0d*asin(1.)    ; Pi\r
+  ZRDSDG= XPI/180.0d          ; Radian to Degree conversion factor\r
+  ZXBM0 = 0.0d\r
+  ZYBM0 = 0.0d\r
+\r
+;=================================================;\r
+  if (ZRPK.eq.0) then\r
+; MERCATOR\r
+;=================================================;\r
+    XBETA=0.\r
+    XLAT0=0.   ; map reference latitude  (degrees)\r
+    ZXBM0 = 0.\r
+    ZYBM0 = 0.\r
+    ZCGAM    = cos(-ZRDSDG*XBETA)\r
+    ZSGAM    = sin(-ZRDSDG*XBETA)\r
+    ZRACLAT0 = XRADIUS*cos(ZRDSDG*ZLAT0)\r
+    do ji=0,IMAX-1\r
+      jj=0\r
+      ZXMI0 = XHAT(ji)-ZXBM0\r
+      ZYMI0 = YHAT(jj)-ZYBM0\r
+      zlon  = (ZXMI0*ZCGAM+ZYMI0*ZSGAM)/(ZRACLAT0*ZRDSDG)+ZLONOR\r
+      do jj=0,JMAX-1\r
+        plon(jj,ji)=zlon\r
+      end do\r
+    end do\r
+    do jj=0,JMAX-1\r
+      ji=0\r
+      ZXMI0 = XHAT(ji)-ZXBM0\r
+      ZYMI0 = YHAT(jj)-ZYBM0\r
+      ZT1  = log(tan(XPI/4.+ZLATOR*ZRDSDG/2.))\r
+      ZT2  = (-ZXMI0*ZSGAM+ZYMI0*ZCGAM)/ZRACLAT0\r
+      zlat  = (-XPI/2.+2.*atan(exp(ZT1+ZT2)))/ZRDSDG\r
+      do ji=0,IMAX-1\r
+        plat(jj,ji)=zlat\r
+      end do\r
+    end do\r
+    \r
+;=================================================;\r
+  else\r
+; STEREOGRAPHIC PROJECTION\r
+;=================================================;\r
+    ZCLAT0  = cos(ZRDSDG*ZLAT0)\r
+    ZSLAT0  = sin(ZRDSDG*ZLAT0)\r
+    ZCLATOR = cos(ZRDSDG*ZLATOR)\r
+    ZSLATOR = sin(ZRDSDG*ZLATOR)\r
+    ZRO0    = (XRADIUS/ZRPK)*(abs(ZCLAT0))^(1.-ZRPK) * \\r
+              ((1.+ZSLAT0)*abs(ZCLATOR)/(1.+ZSLATOR))^ZRPK\r
+    ZGA0    = (ZRPK*(ZLONOR-ZLON0)-ZBETA)*ZRDSDG\r
+    ZXP     = ZXBM0-ZRO0*sin(ZGA0)\r
+    ZYP     = ZYBM0+ZRO0*cos(ZGA0)\r
+    do ji=0,IMAX-1\r
+      do jj=0,JMAX-1\r
+        ZATA = atan2( -(ZXP-XHAT(ji)) , (ZYP-YHAT(jj)) )/ZRDSDG\r
+        zlon  = (ZBETA+ZATA)/ZRPK+ZLON0\r
+        plon(jj,ji)=zlon\r
+        ZRO2 = (XHAT(ji)-ZXP)^2+(YHAT(jj)-ZYP)^2\r
+        ZJD1 = XRADIUS*(abs(ZCLAT0))^(1.-ZRPK)\r
+        ZT1  = (ZJD1)^(2./ZRPK)* (1+ZSLAT0)^2\r
+        ZJD3 = (ZRPK^2*ZRO2)\r
+        ZT2  = ZJD3\r
+        ZT2 = ZT2^(1./ZRPK)\r
+        ZJD1 = (ZT1-ZT2)/(ZT1+ZT2)\r
+        ZJD1 = acos(ZJD1)\r
+        ZJD3 = ZJD1\r
+        zlat = (XPI/2.-ZJD3)/ZRDSDG\r
+        plat(jj,ji)=zlat\r
+      end do\r
+    end do\r
+    \r
+  end if\r
+\r
+; Defining the corners of the domain\r
+;====================================\r
+  if (icorner(0,0).eq.icorner(1,1)) then\r
+    icorner(0,0)=0\r
+    icorner(1,0)=JMAX-1\r
+    icorner(0,1)=0\r
+    icorner(1,1)=IMAX-1\r
+  end if\r
+; print ("icorner"+icorner)\r
+\r
+  res@mpLimitMode            = "Corners"\r
+  res@mpLeftCornerLatF       = plat(icorner(0,0),icorner(0,1))\r
+  res@mpLeftCornerLonF       = plon(icorner(0,0),icorner(0,1))\r
+  res@mpRightCornerLatF     = plat(icorner(1,0),icorner(1,1))\r
+  res@mpRightCornerLonF     = plon(icorner(1,0),icorner(1,1))\r
+\r
+; print ("Corner (0,0); Lat="+res@mpLeftCornerLatF+ \\r
+;                    ", Lon="+res@mpLeftCornerLonF)\r
+; print ("Oppos corner; Lat="+res@mpRightCornerLatF+ \\r
+;                     ", Lon= "+res@mpRightCornerLonF)\r
+\r
+;==========================================\r
+; Turn on lat / lon labeling\r
+;==========================================\r
+  res@pmTickMarkDisplayMode = "Always"         ; turn on tickmarks\r
+  res@mpOutlineBoundarySets = "AllBoundaries" ; state boundaries\r
+  res@mpPerimDrawOrder      = "PostDraw"       ; force map perim\r
+;==========================================\r
+; Needed for regional native projection \r
+;==========================================\r
+  res@tfDoNDCOverlay       = True      \r
+  res@gsnAddCyclic         = False              ; regional data\r
+  \r
+end\r
+\r
+;===========================================\r
+;------------------------------------------------------------------------\r
+undef("MESONH_pinter")\r
+function MESONH_pinter( pfield:numeric, loc_param:numeric, ppabs:numeric )\r
+;*************************************************************************\r
+; S. BIELLI\r
+; This is a routine that interpolate  fields on pressure level for plotting\r
+; based on pinter.f90 \r
+; The field to be interpolated must be given at the mass point (grid 1)\r
+; usage : var_inter=MESONHfunction(var_to_interpol, 850., AbsPressure)\r
+; Abs pressure must be in Pa\r
+;\r
+\r
+begin\r
+\r
+ dimL= dimsizes(loc_param)\r
+\r
+; First test for grid = 0\r
+\r
+  dimp=dimsizes(ppabs)\r
+\r
+  pout=pfield(0:dimL-1,:,:)\r
+  pfield@_FillValue=999\r
+  pout@_FillValue=999\r
+  pout=pout@_FillValue\r
+\r
+  do jkp = 0, dimL-1\r
+       zref=log10(loc_param(jkp)*100.)\r
+       do jloop = 0, dimp(1)-1\r
+         do iloop = 0, dimp(2)-1\r
+            kloop=0\r
+            flag=True\r
+           do while (flag .and. (kloop.lt.(dimp(2)-2)))\r
+             if (.not.ismissing(ppabs(kloop,jloop,iloop))) then\r
+               zxm=log10(ppabs(kloop,jloop,iloop))\r
+               zxp=log10(ppabs(kloop+1,jloop,iloop))\r
+               if ((zxp-zref)*(zref-zxm) .ge. 0) then\r
+                  pout(jkp,jloop,iloop)= (pfield(kloop,jloop,iloop)*(zxp-zref)+ \\r
+                       pfield(kloop+1,jloop,iloop)*(zref-zxm))/ (zxp-zxm)\r
+                  flag=False\r
+                end if\r
+              end if\r
+              kloop=kloop+1\r
+            end do\r
+          end do\r
+        end do\r
+  end do\r
+\r
+  return(pout)\r
+       \r
+end\r
+\r
+;--------------------------------------------------------------------------------\r
+undef("mnh_map")\r
+function mnh_map(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical)   \r
+\r
+begin\r
+;\r
+; This function creates a map plot, and bases the projection on\r
+; the MAP_PROJ attribute in the given file.\r
+;\r
+;   1. Make a copy of the resource list, and set some resources\r
+;      common to all map projections.\r
+;\r
+;   2. Determine the projection being used, and set resources based\r
+;      on that projection.\r
+;\r
+;   3. Create the map plot, and draw and advance the frame\r
+;      (if requested).\r
+\r
+    opts = opt_args      ; Make a copy of the resource list\r
+    opts  = True\r
+\r
+; Set some resources depending on what kind of map projection is \r
+; chosen.\r
+;\r
+;   ZRPK != 0 : "Stereographic"\r
+;   ZRPK = 0 : "Mercator"\r
+;=================================================;\r
+; src/mesonh_MOD/mode_gridproj.f90\r
+;=================================================;\r
+  XRADIUS=6371229.0d ; Earth radius (meters)\r
+  XPI=2.0d*asin(1.)    ; Pi\r
+  ZRDSDG= XPI/180.0d          ; Radian to Degree conversion factor\r
+  ZXBM0 = 0.0d\r
+  ZYBM0 = 0.0d\r
+\r
+  if(isfilevar(in_file,"RPK"))\r
+    ZRPK=in_file->RPK\r
+    ZLON0=in_file->LON0\r
+    ZLAT0=in_file->LAT0\r
+    ZLATOR=in_file->LATOR\r
+    ZLONOR=in_file->LONOR\r
+    ZBETA=in_file->BETA\r
+  else\r
+   print ("mnh_map: Error no RPK variable in input file")\r
+   return(new(1,graphic))\r
+  end if\r
\r
+; Case netcdf from lfi2cdf\r
+   if(isfilevar(in_file,"IMAX"))\r
+       XHAT=in_file->XHAT\r
+       YHAT=in_file->YHAT\r
+       IMAX= dimsizes(XHAT)-2\r
+       JMAX= dimsizes(YHAT)-2\r
+       zdx=XHAT(2)-XHAT(1)\r
+       zdy=YHAT(2)-YHAT(1)\r
+       do ji=0,IMAX-1\r
+         XHAT(ji)=XHAT(ji)+zdx*1.5\r
+       end do\r
+       do jj=0,JMAX-1\r
+         YHAT(jj)=YHAT(jj)+zdy*1.5\r
+       end do\r
+   else\r
+; Case netcdf from extractdia\r
+       XHAT=in_file->W_E_direction\r
+       YHAT=in_file->S_N_direction\r
+       IMAX= dimsizes(XHAT)\r
+       JMAX= dimsizes(YHAT)\r
+   end if\r
+;\r
+  \r
+  lat = new((/JMAX,IMAX/),"double")\r
+  lon = new((/JMAX,IMAX/),"double")\r
+\r
+\r
+;   Stereographic projection\r
+      if(ZRPK .gt. 0)\r
+        projection          = "Stereographic"\r
+        opts@mpCenterLatF  = get_res_value_keep(opts, "mpCenterLatF", 90)\r
+        opts@mpCenterLonF  = get_res_value_keep(opts, "mpCenterLonF",ZLON0)\r
+        opts@mpCenterRotF  = get_res_value_keep(opts, "mpCenterRotF",ZBETA)\r
+      end if\r
+\r
+      if(ZRPK .lt. 0)\r
+        projection          = "Stereographic"\r
+        opts@mpCenterLatF  = get_res_value_keep(opts, "mpCenterLatF", -90)\r
+        opts@mpCenterLonF  = get_res_value_keep(opts, "mpCenterLonF",ZLON0)\r
+        opts@mpCenterRotF  = get_res_value_keep(opts, "mpCenterRotF",ZBETA)\r
+      end if\r
+\r
+;   Mercator projection\r
+      if(ZRPK .eq. 0)\r
+        projection          = "Mercator"\r
+      end if\r
+\r
+    opts@mpNestTime = get_res_value_keep(opts, "mpNestTime",0)\r
+\r
+    \r
+; LAT and LON are not saved in the file \r
+  if (ZRPK.eq.0) then\r
+    XBETA=0.\r
+    XLAT0=0.   ; map reference latitude  (degrees)\r
+    ZXBM0 = 0.\r
+    ZYBM0 = 0.\r
+    ZCGAM    = cos(-ZRDSDG*XBETA)\r
+    ZSGAM    = sin(-ZRDSDG*XBETA)\r
+    ZRACLAT0 = XRADIUS*cos(ZRDSDG*ZLAT0)\r
+    do ji=0,IMAX-1\r
+      jj=0\r
+      ZXMI0 = XHAT(ji)-ZXBM0\r
+      ZYMI0 = YHAT(jj)-ZYBM0\r
+      zlon  = (ZXMI0*ZCGAM+ZYMI0*ZSGAM)/(ZRACLAT0*ZRDSDG)+ZLONOR\r
+      do jj=0,JMAX-1\r
+        lon(jj,ji)=zlon\r
+      end do\r
+    end do\r
+    do jj=0,JMAX-1\r
+      ji=0\r
+      ZXMI0 = XHAT(ji)-ZXBM0\r
+      ZYMI0 = YHAT(jj)-ZYBM0\r
+      ZT1  = log(tan(XPI/4.+ZLATOR*ZRDSDG/2.))\r
+      ZT2  = (-ZXMI0*ZSGAM+ZYMI0*ZCGAM)/ZRACLAT0\r
+      zlat  = (-XPI/2.+2.*atan(exp(ZT1+ZT2)))/ZRDSDG\r
+      do ji=0,IMAX-1\r
+        lat(jj,ji)=zlat\r
+      end do\r
+    end do\r
+   else\r
+    ZCLAT0  = cos(ZRDSDG*ZLAT0)\r
+    ZSLAT0  = sin(ZRDSDG*ZLAT0)\r
+    ZCLATOR = cos(ZRDSDG*ZLATOR)\r
+    ZSLATOR = sin(ZRDSDG*ZLATOR)\r
+    ZRO0    = (XRADIUS/ZRPK)*(abs(ZCLAT0))^(1.-ZRPK) * \\r
+              ((1.+ZSLAT0)*abs(ZCLATOR)/(1.+ZSLATOR))^ZRPK\r
+    ZGA0    = (ZRPK*(ZLONOR-ZLON0)-ZBETA)*ZRDSDG\r
+    ZXP     = ZXBM0-ZRO0*sin(ZGA0)\r
+    ZYP     = ZYBM0+ZRO0*cos(ZGA0)\r
+    do ji=0,IMAX-1\r
+      do jj=0,JMAX-1\r
+        ZATA = atan2( -(ZXP-XHAT(ji)) , (ZYP-YHAT(jj)) )/ZRDSDG\r
+        zlon  = (ZBETA+ZATA)/ZRPK+ZLON0\r
+        lon(jj,ji)=zlon\r
+        ZRO2 = (XHAT(ji)-ZXP)^2+(YHAT(jj)-ZYP)^2\r
+        ZJD1 = XRADIUS*(abs(ZCLAT0))^(1.-ZRPK)\r
+        ZT1  = (ZJD1)^(2./ZRPK)* (1+ZSLAT0)^2\r
+        ZJD3 = (ZRPK^2*ZRO2)\r
+        ZT2  = ZJD3\r
+        ZT2 = ZT2^(1./ZRPK)\r
+        ZJD1 = (ZT1-ZT2)/(ZT1+ZT2)\r
+        ZJD1 = acos(ZJD1)\r
+        ZJD3 = ZJD1\r
+        zlat = (XPI/2.-ZJD3)/ZRDSDG\r
+        lat(jj,ji)=zlat\r
+      end do\r
+    end do\r
+    end if\r
+\r
+      dims = dimsizes(lat)\r
+\r
+      do ii = 0, dims(0)-1\r
+      do jj = 0, dims(1)-1\r
+        if ( lon(ii,jj) .lt. 0.0) then\r
+          lon(ii,jj) = lon(ii,jj) + 360.\r
+        end if\r
+      end do\r
+      end do\r
+\r
+      opts@start_lat = lat(0,0)\r
+      opts@start_lon = lon(0,0)\r
+      opts@end_lat   = lat(dims(0)-1,dims(1)-1)\r
+      opts@end_lon   = lon(dims(0)-1,dims(1)-1)\r
+\r
+\r
+; Set some resources common to all map projections.\r
+  opts = set_mp_resources(opts)\r
+\r
+      if ( isatt(opts,"ZoomIn") .and. opts@ZoomIn ) then\r
+        y1 = 0\r
+        x1 = 0\r
+        y2 = dims(0)-1\r
+        x2 = dims(1)-1\r
+        if ( isatt(opts,"Ystart") ) then\r
+          y1 = opts@Ystart\r
+          delete(opts@Ystart)\r
+        end if\r
+        if ( isatt(opts,"Xstart") ) then\r
+          x1 = opts@Xstart\r
+          delete(opts@Xstart)\r
+        end if\r
+        if ( isatt(opts,"Yend") ) then\r
+          if ( opts@Yend .le. y2 ) then\r
+            y2 = opts@Yend\r
+          end if\r
+          delete(opts@Yend)\r
+        end if\r
+        if ( isatt(opts,"Xend") ) then\r
+          if ( opts@Xend .le. x2 ) then\r
+            x2 = opts@Xend\r
+          end if\r
+          delete(opts@Xend)\r
+        end if\r
+\r
+        opts@mpLeftCornerLatF      = lat(y1,x1)\r
+        opts@mpLeftCornerLonF      = lon(y1,x1)\r
+        opts@mpRightCornerLatF     = lat(y2,x2)\r
+        opts@mpRightCornerLonF     = lon(y2,x2)\r
+\r
+        if ( opts@mpRightCornerLonF .lt. 0.0 ) then\r
+          opts@mpRightCornerLonF  = opts@mpRightCornerLonF + 360.0\r
+        end if \r
+\r
+        delete(opts@ZoomIn)\r
+      end if\r
+\r
+\r
+; The default is not to draw the plot or advance the frame, and\r
+; to maximize the plot in the frame.\r
+\r
+  opts@gsnDraw       = get_res_value_keep(opts,"gsnDraw",     False)\r
+  opts@gsnFrame      = get_res_value_keep(opts,"gsnFrame",    False)\r
+  opts@gsnMaximize   = get_res_value_keep(opts,"gsnMaximize", True)\r
+\r
+  delete_attrs(opts)                             ; Clean up.\r
+  mp = gsn_map(wks,projection,opts)              ; Create map plot.\r
+\r
+  return(mp)                                     ; Return.\r
+\r
+end\r
+\r
+;--------------------------------------------------------------------------------\r
+\r
+undef("mnh_map_overlays")\r
+function mnh_map_overlays(in_file[1]:file, \\r
+                          wks:graphic, \\r
+                          plots[*]:graphic, \\r
+                          opt_arg[1]:logical, \\r
+                          opt_mp[1]:logical) \r
+\r
+; Based on wrf_map_overlays \r
+; \r
+; This procedure takes an array of plots and overlays them on a\r
+; base plot - map background.\r
+;\r
+; It will advance the plot and cleanup, unless you set the\r
+; PanelPlot resource to True.\r
+;\r
+;  Attributes recognized by this procedure:\r
+;     FramePlot\r
+;     PanelPlot\r
+;     NoTitles                  (don't do any titles) \r
+;     CommonTitle & PlotTile is used to overwrite field titles\r
+;        CommonTitle will super-seed NoTitles\r
+;\r
+; If FramePlot False, then Draw the plot but do not Frame.\r
+; In this case a user want to add to the drawing, and will\r
+; have to advance the Frame manually in the script.\r
+;\r
+; If the "NoTitles" attribute exists and is set True, then\r
+; don't create the top-left titles, and leave the main titles alone.\r
+; This resource can be useful if you are planning to panel\r
+; the plots.\r
+;\r
+; If PanelPlot is set to True, then this flags to wrf_map_overlays\r
+; that these plots are going to be eventually paneled (likely\r
+; by gsn_panel), and hence 1) draw and frame should not be called\r
+; (unless gsnDraw and/or gsnFrame are explicitly set to True),\r
+; and 2) the overlays and titles should not be removed with\r
+; NhlRemoveOverlay and NhlRemoveAnnotation.\r
+;\r
+begin\r
+\r
+  opts     = opt_arg         ; Make a copy of the resource lists\r
+  opt_mp_2 = opt_mp      \r
+\r
+  ; Let's make the map first\r
+  base = mnh_map(wks,in_file,opt_mp_2)\r
+\r
+  no_titles  = get_res_value(opts,"NoTitles",False)     ; Do we want field titles?\r
+  com_title  = get_res_value(opts,"CommonTitle",False)     ; Do we have a common title?\r
+  if ( com_title ) then\r
+    plot_title = get_res_value(opts,"PlotTitle","  ")\r
+    no_titles = True\r
+  end if\r
+  \r
+  call_draw  = True\r
+  call_frame = get_res_value(opts,"FramePlot",True)     ; Do we want to frame the plot?\r
+  panel_plot = get_res_value(opts,"PanelPlot",False)    ; Are we paneling?\r
+  opts@gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)\r
+\r
+  nplots = dimsizes(plots)\r
+;  font_color = "Black"\r
+\r
+  do i=0,nplots-1\r
+    if(.not.ismissing(plots(i))) then\r
+;      class_name = NhlClassName(plots(i))\r
+;      print(class_name)\r
+;      if(class_name.eq."contourPlotClass") then\r
+;        getvalues plots(i)\r
+;          "cnFillOn"    : fill_on\r
+;          "cnLineColor" : line_color\r
+;        end getvalues\r
+;        if (.not.fill_on) then\r
+;          font_color = line_color      \r
+;        end if       \r
+;      end if\r
+      if(.not.no_titles) then\r
+        getvalues plots(i)\r
+          "tiMainString" : SubTitle\r
+        end getvalues\r
+        if(i.eq.0) then\r
+          SubTitles = SubTitle\r
+        else\r
+          SubTitles = SubTitles + "~C~" + SubTitle\r
+        end if\r
+      end if\r
+      if(com_title .and. i .eq. nplots-1) then\r
+        getvalues plots(i)\r
+          "tiMainString" : SubTitle\r
+        end getvalues\r
+        SubTitles = plot_title\r
+      end if\r
+      setvalues plots(i)\r
+        "tfDoNDCOverlay" : True\r
+        "tiMainOn"       : False\r
+      end setvalues\r
+      overlay(base,plots(i))\r
+    else\r
+      print("mnh_map_overlays: Warning: overlay plot #" + i + " is not valid.")\r
+    end if\r
+  end do\r
+\r
+  if(.not.no_titles .or. com_title) then\r
+    font_height = get_res_value_keep(opts,"FontHeightF",0.01)\r
+    txt = create "map_titles" textItemClass wks\r
+      "txString"      : SubTitles\r
+      "txFontHeightF" : font_height\r
+     ;"txFontColor"   : font_color\r
+    end create\r
+    anno = NhlAddAnnotation(base,txt)\r
+    setvalues anno\r
+      "amZone"           : 3\r
+      "amJust"           : "BottomLeft"\r
+      "amSide"           : "Top"\r
+      "amParallelPosF"   : 0.005\r
+      "amOrthogonalPosF" : 0.03\r
+      "amResizeNotify"   : False\r
+    end setvalues\r
+    base@map_titles = anno\r
+  end if\r
+;\r
+; gsnDraw and gsnFrame default to False if panel plot.\r
+;\r
+  if(panel_plot) then\r
+    call_draw = False\r
+    call_frame= False\r
+  end if\r
+\r
+\r
+  opts@gsnDraw     = get_res_value_keep(opts,"gsnDraw",     call_draw)\r
+  opts@gsnFrame    = get_res_value_keep(opts,"gsnFrame",    call_frame)\r
+\r
+  draw_and_frame(wks,base,opts@gsnDraw,opts@gsnFrame,False, \\r
+                 opts@gsnMaximize)\r
+\r
+  if(.not.panel_plot) then\r
+    do i=0,nplots-1\r
+      if(.not.ismissing(plots(i))) then\r
+        NhlRemoveOverlay(base,plots(i),False)\r
+      else\r
+        print("wrf_remove_map_overlays: Warning: overlay plot #" + i + " is not valid.")\r
+        print("                        Nothing to remove.")\r
+      end if\r
+    end do\r
+  end if\r
+\r
+  if(.not.no_titles.and..not.panel_plot) then\r
+    if(isatt(base,"map_titles")) then\r
+      NhlRemoveAnnotation(base,base@map_titles)\r
+      delete(base@map_titles)\r
+    end if\r
+  end if\r
+\r
+return(base)\r
+end\r
+\r
+;--------------------------------------------------------------------------------\r
+undef("wrf_user_intrp3d")\r
+function wrf_user_intrp3d( var3d:numeric, z_in:numeric, \\r
+                           plot_type:string, \\r
+                           loc_param:numeric, angle:numeric, opts:logical )\r
+\r
+; var3d      - 3d field to interpolate (all input fields must be unstaggered)\r
+; z_in       - interpolate to this field (either p/z)\r
+; plot_type  - interpolate horizontally "h", or vertically "v"\r
+; loc_param  - level(s) for horizontal plots (eg. 500hPa ; 3000m - scalar),\r
+;              plane for vertical plots (2 values representing an xy point\r
+;              on the model domain through which the vertical plane will pass\r
+;              OR 4 values specifying start and end values\r
+; angle      - 0.0 for horizontal plots, and\r
+;              an angle for vertical plots - 90 represent a WE cross section\r
+; opts         Used IF opts is TRUE, else use loc_param and angle to determine crosssection\r
+\r
+begin\r
+\r
+\r
+  if(plot_type .eq. "h" ) then   ;  horizontal cross section needed\r
+\r
+     dimL = dimsizes(loc_param)\r
+\r
+     dims = dimsizes(var3d)\r
+     nd = dimsizes(dims)\r
+\r
+     dimX = dims(nd-1)\r
+     dimY = dims(nd-2)\r
+     dimZ = dims(nd-3)\r
+     dim4 = 1\r
+     dim5 = 1\r
+     if ( nd .eq. 4 ) then\r
+       dim4 = dims(nd-4)\r
+     end if\r
+     if ( nd .eq. 5 ) then\r
+       dim4 = dims(nd-4)\r
+       dim5 = dims(nd-5)\r
+     end if\r
+\r
+     var3  = new ( (/ dim5, dim4, dimZ, dimY, dimX /) , typeof(var3d) )\r
+     z     = new ( (/ dim5, dim4, dimZ, dimY, dimX /) , typeof(var3d) )\r
+     var2d = new ( (/ dim5, dim4, dimL, dimY, dimX /) , typeof(var3d) )\r
+\r
+     if ( nd .eq. 5 ) then\r
+       var3 = var3d\r
+          z =  z_in\r
+     end if\r
+     if ( nd .eq. 4 ) then\r
+       var3(0,:,:,:,:) = var3d(:,:,:,:)\r
+          z(0,:,:,:,:) =  z_in(:,:,:,:)\r
+     end if\r
+     if ( nd .eq. 3 ) then\r
+       var3(0,0,:,:,:) = var3d(:,:,:)\r
+          z(0,0,:,:,:) =  z_in(:,:,:)\r
+     end if\r
+\r
+\r
+     if ( z(0,0,0,0,0) .gt. 500.) then\r
+       ; We must be interpolating to pressure\r
+       ; This routine needs input field and level in hPa - lets make sure of this\r
+       if ( z(0,0,0,0,0) .gt. 2000. ) then\r
+         ; looks like we have Pa as input - make this hPa\r
+         z = z * 0.01\r
+       end if\r
+       if ( loc_param(0) .gt. 2000. ) then\r
+         ; looks like the input was specified in Pa - change this\r
+         loc_param = loc_param * 0.01\r
+       end if\r
+     end if\r
+\r
+     do il = 0,dimL-1\r
+       var = wrf_interp_3d_z(var3,z,loc_param(il))\r
+       var2d(:,:,il,:,:) = var(:,:,:,:)\r
+     end do\r
+\r
+     copy_VarAtts(var3d,var3)\r
+     if(isatt(var3,"description")) then\r
+       delete_VarAtts(var3,(/"description"/))\r
+     end if\r
+     if(isatt(var3,"units")) then\r
+       delete_VarAtts(var3,(/"units"/))\r
+     end if\r
+     if(isatt(var3,"MemoryOrder")) then\r
+       delete_VarAtts(var3,(/"MemoryOrder"/))\r
+     end if\r
+     if(isatt(var3,"_FillValue")) then\r
+       delete_VarAtts(var3,(/"_FillValue"/))\r
+     end if\r
+     copy_VarAtts(var3,var2d)\r
+\r
+     nn = nd-2\r
+     var2d!nn = "plevs"\r
+\r
+     if ( dimL .gt. 1 ) then\r
+       if ( nd .eq. 5 ) then\r
+         return( var2d )\r
+       end if\r
+       if ( nd .eq. 4 ) then\r
+         return( var2d(0,:,:,:,:) )\r
+       end if\r
+       if ( nd .eq. 3 ) then\r
+         return( var2d(0,0,:,:,:) )\r
+       end if\r
+     else\r
+       if ( z(0,0,0,0,0) .gt. 500.) then\r
+          var2d@PlotLevelID = loc_param + " hPa"\r
+       else\r
+          var2d@PlotLevelID = .001*loc_param + " km"\r
+       end if\r
+       if ( nd .eq. 5 ) then\r
+         return( var2d(:,:,0,:,:) )\r
+       end if\r
+       if ( nd .eq. 4 ) then\r
+         return( var2d(0,:,0,:,:) )\r
+       end if\r
+       if ( nd .eq. 3 ) then\r
+         return( var2d(0,0,0,:,:) )\r
+       end if\r
+     end if\r
+\r
+\r
+  end if\r
+\r
+\r
+\r
+\r
+  if(plot_type .eq. "v" ) then   ;  vertical cross section needed\r
+\r
+     dims = dimsizes(var3d)\r
+     if ( dimsizes(dims) .eq. 4 ) then\r
+       if ( z_in(0,0,0,0) .gt. 500.) then\r
+         ; We must be interpolating to pressure\r
+           ; This routine needs input field and level in hPa - lets make sure of this\r
+         if ( z_in(0,0,0,0) .gt. 2000. ) then\r
+           ; looks like we have Pa as input - make this hPa\r
+           z_in = z_in * 0.01\r
+         end if\r
+       end if\r
+       z = z_in(0,:,:,:)\r
+     else\r
+       if ( z_in(0,0,0) .gt. 500.) then\r
+         ; We must be interpolating to pressure\r
+           ; This routine needs input field and level in hPa - lets make sure of this\r
+         if ( z_in(0,0,0) .gt. 2000. ) then\r
+           ; looks like we have Pa as input - make this hPa\r
+           z_in = z_in * 0.01\r
+         end if\r
+       end if\r
+       z = z_in\r
+     end if\r
+\r
+; set vertical cross section\r
+     if (opts) then\r
+       xy = wrf_user_set_xy( z, loc_param(0)-1, loc_param(1)-1, \    ; the -1 is for NCL dimensions\r
+                                loc_param(2)-1, loc_param(3)-1, \\r
+                                angle, opts )\r
+     else\r
+       xy = wrf_user_set_xy( z, loc_param(0), loc_param(1), \\r
+                                0.0, 0.0, angle, opts )\r
+     end if\r
+     xp = dimsizes(xy)\r
+\r
+\r
+; first we interp z\r
+     var2dz   = wrf_interp_2d_xy( z, xy)\r
+\r
+;  interp to constant z grid\r
+     if(var2dz(0,0) .gt. var2dz(1,0) ) then  ; monotonically decreasing coordinate\r
+        z_max = floor(max(z)/10)*10     ; bottom value\r
+        z_min = ceil(min(z)/10)*10      ; top value\r
+        dz = 1.\r
+        nlevels = tointeger( (z_max-z_min)/dz)\r
+        z_var2d = new( (/nlevels/), typeof(z))\r
+        z_var2d(0) = z_max\r
+        dz = -dz\r
+     else\r
+        z_max = max(z)\r
+        z_min = 0.\r
+;; MODI SOLINE\r
+;        dz = 0.01 * z_max\r
+        dz = 0.001 * z_max\r
+        nlevels = tointeger( z_max/dz )\r
+        z_var2d = new( (/nlevels/), typeof(z))\r
+        z_var2d(0) = z_min\r
+     end if\r
+;     print("nlevels="+nlevels)\r
+;     print("dz="+dz)\r
+\r
+     do i=1, nlevels-1\r
+        z_var2d(i) = z_var2d(0)+i*dz\r
+     end do\r
+\r
+\r
+; interp the variable\r
+     if ( dimsizes(dims) .eq. 4 ) then\r
+       var2d = new( (/dims(0), nlevels, xp(0)/), typeof(var2dz))\r
+       do it = 0,dims(0)-1\r
+         var2dtmp = wrf_interp_2d_xy( var3d(it,:,:,:), xy)\r
+         do i=0,xp(0)-1\r
+            var2d(it,:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d)\r
+         end do\r
+       end do\r
+       var2d!0 = var3d!0\r
+       var2d!1 = "Vertical"\r
+       var2d!2 = "Horizontal"\r
+     else\r
+       var2d = new( (/nlevels, xp(0)/), typeof(var2dz))\r
+       var2dtmp = wrf_interp_2d_xy( var3d, xy)\r
+       do i=0,xp(0)-1\r
+          var2d(:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d)\r
+       end do\r
+       var2d!0 = "Vertical"\r
+       var2d!1 = "Horizontal"\r
+     end if\r
+\r
+\r
+     st_x = tointeger(xy(0,0)) + 1\r
+     st_y = tointeger(xy(0,1)) + 1\r
+     ed_x = tointeger(xy(xp(0)-1,0)) + 1\r
+     ed_y = tointeger(xy(xp(0)-1,1)) + 1\r
+     if (opts) then\r
+       var2d@Orientation = "Cross-Sesion: (" + \\r
+                            st_x + "," + st_y + ") to (" + \\r
+                            ed_x + "," + ed_y + ")"\r
+     else\r
+       var2d@Orientation = "Cross-Sesion: (" + \\r
+                            st_x + "," + st_y + ") to (" + \\r
+                            ed_x + "," + ed_y + ") ; center=(" + \\r
+                            loc_param(0) + "," + loc_param(1) + \\r
+                            ") ; angle=" + angle\r
+     end if\r
+\r
+     return(var2d)\r
+end if\r
+\r
+\r
+end\r
+\r
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/clean_ncl b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/clean_ncl
new file mode 100755 (executable)
index 0000000..9fb90e5
--- /dev/null
@@ -0,0 +1,16 @@
+echo "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
+echo "!!!!  WARNING :: CLEAN_RANGS_FILES       !!!!"
+echo "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
+echo
+echo do you real want to delete high resolution
+echo
+echo coastline files for ncl  ?????
+echo
+echo "(yes/no)"
+read ANSWER
+set -x
+if [ "X$ANSWER" = "Xyes" ]
+then
+rm -fr  ran* gs*
+fi
+rm -f *.lfi *.nc4 
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/get_ncl_highres_files b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/get_ncl_highres_files
new file mode 100755 (executable)
index 0000000..b780817
--- /dev/null
@@ -0,0 +1,20 @@
+NCARG_RANGS=`pwd`
+export NCARG_RANGS
+echo 
+[ -f ~/.hluresfile ] || ( wget -c -nd http://www.ncl.ucar.edu/Document/Graphics/hluresfile ; mv hluresfile ~/.hluresfile )
+echo "You will download high resolution coastline files for NCL"
+echo "in this directory:"
+echo NCL_HIGHRES_FILES=$NCARG_RANGS 
+echo 
+set -x
+cd $NCARG_RANGS
+NCL_HIGHRES_URL="http://www2008.io-warnemuende.de/homepages/rfeistel/download"
+export NCL_HIGHRES_URL
+for file in rangs\(0\) rangs\(1\) rangs\(2\) \
+            rangs\(3\) rangs\(4\)  \
+            gshhs\(0\) gshhs\(1\) gshhs\(2\) \
+            gshhs\(3\) gshhs\(4\)
+do
+[ -f $file.zip ] || ( wget -c -nd $NCL_HIGHRES_URL/$file.zip ; unzip $file.zip ; ) 
+done
+
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_BasicMap.ncl b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_BasicMap.ncl
new file mode 100644 (file)
index 0000000..c8f4b66
--- /dev/null
@@ -0,0 +1,179 @@
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
+
+;=============================================================;
+;  Example of script to plot ZS field from MESONH netcdf4 file
+;  Could be use with any files generated by MESONH if LCARTESIAN=False
+;=============================================================;
+begin
+;=============================================================;
+; Open file 
+;=============================================================;
+
+  mnh_file="REUNI.1.00A20.004.nc4"
+  a = addfile(mnh_file, "r")
+
+;=================================================;
+; Get information on variable sizes
+;=================================================;
+
+  mdims = getfilevardimsizes(a,"ZS") ; get some dimension sizes for the file
+  nd = dimsizes(mdims)
+  imax=mdims(nd-1)-2  ; -2 : to remove non-physical values
+  jmax=mdims(nd-2)-2  ; -2 : to remove non-physical values
+
+;=================================================;
+; Read the variables 
+;=================================================;
+
+; Note: do not read first and last value which are non physical values
+; -------------------------------------
+
+  lat2d=a->LAT(1:jmax,1:imax)
+  lon2d=a->LON(1:jmax,1:imax)
+
+  ZS = a->ZS(1:jmax,1:imax)
+  ZS@long_name="Orography"
+  ZS@units="m"
+
+; Associate coordinates to variable
+; Necessary to plot the field at the right geographical position 
+; -------------------------------------
+
+  ZS@lat2d = lat2d
+  ZS@lon2d = lon2d
+
+; Read projection parameters 
+; -------------------------------------
+
+  RPK = a->RPK
+  BETA =a->BETA
+  LON0= a->LON0
+  
+;=================================================;
+; CREATE PLOT : Orography (ZS, filled contour)
+;=================================================;
+; Open workstation and define colormap
+; -------------------------------------
+
+  type  = "x11" ; open a x11 window
+                ; change type to png, ps or pdf to get the plot into a file
+  wks  = gsn_open_wks(type,"plt_BasicMap")   ; 
+  gsn_define_colormap(wks,"topo_15lev") ; Choose colormap
+
+
+;==========================================
+; Map ressources
+;==========================================
+  resmap                        = True         
+
+; If there is an error on HighRes, it means that you don't have the HighRes data
+; You need to download them or the change HighRes by MediumRes
+; See https://www.ncl.ucar.edu/Document/Graphics/rangs.shtml for info
+; -------------------------------------
+  resmap@mpDataBaseVersion     = "HighRes"     ; choose highres map data version (must be donwloaded)
+  resmap@mpGridAndLimbOn       = True             ; turn on lat/lon lines
+  resmap@mpGridLatSpacingF     = 10.              ; spacing for lat lines
+  resmap@mpGridLonSpacingF     = 10.              ; spacing for lon lines
+
+  resmap@pmTickMarkDisplayMode = "Always"         ; turn on tickmarks
+
+;=================================================;
+; Set map projection ressources using projection parameters
+;=================================================;
+  if (RPK.gt.0)
+; ---------------------------
+  ;   Lambert  projection from north pole
+; ---------------------------
+   resmap@mpProjection          = "LambertConformal"     ; projection
+   pole                         = 1    ; projection for north hemisphere 
+   resmap@mpLambertParallel1F   = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere
+   resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file
+   resmap@mpLambertMeridianF    = LON0      ; ncl adds from grib file
+  end if
+
+  if (RPK.lt.0)
+; ---------------------------
+  ;   Lambert projection from south pole
+; ---------------------------
+   resmap@mpProjection          = "LambertConformal"     ; projection
+   pole                         = -1                     ; projection for south hemisphere
+   resmap@mpLambertParallel1F   = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere
+   resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file
+   resmap@mpLambertMeridianF    =  LON0      ; ncl adds from grib file
+  end if
+
+  if (RPK.eq.1)
+; ---------------------------
+  ;   Stereographic projection north
+; ---------------------------
+    resmap@mpProjection = "Stereographic"
+    resmap@mpCenterLonF           = LON0
+    resmap@mpCenterRotF           = BETA
+    resmap@mpCenterLatF           = 90
+  end if
+
+  if (RPK.eq.-1)
+; ---------------------------
+  ;   Stereographic projection south 
+; ---------------------------
+    resmap@mpProjection = "Stereographic"
+    resmap@mpCenterLonF           = LON0
+    resmap@mpCenterRotF           = BETA
+    resmap@mpCenterLatF           = -90
+  end if
+
+  if (RPK.eq.0) then
+; ---------------------------
+  ;   Mercator projection
+; ---------------------------
+    resmap@mpProjection = "Mercator"
+  end if
+
+ print("Map projection="+resmap@mpProjection)
+
+;==========================================
+; Defining the corners for projection otherwise the plot will be global
+;==========================================
+  resmap@mpLimitMode           = "Corners"
+  resmap@mpLeftCornerLatF      = lat2d(0,0)
+  resmap@mpLeftCornerLonF      = lon2d(0,0)
+  resmap@mpRightCornerLatF     = lat2d(jmax-1,imax-1)
+  resmap@mpRightCornerLonF     = lon2d(jmax-1,imax-1)
+
+
+ print ("Corner (0,0); Lat="+resmap@mpLeftCornerLatF+ \
+                    ", Lon="+resmap@mpLeftCornerLonF)
+ print ("Oppos corner; Lat="+resmap@mpRightCornerLatF+ \
+                     ", Lon= "+resmap@mpRightCornerLonF)
+
+
+;==========================================
+; Create ZS plot  (contour)
+;==========================================
+; ---------------------------
+; General ressources
+; ---------------------------
+  resmap@gsnMaximize            = True          ; maximize size
+  resmap@gsnSpreadColors       = True          ; use full range of colormap
+
+; ---------------------------
+; Contour ressources
+; ---------------------------
+  resmap@cnFillOn              = True          ; turn on color fill
+  resmap@cnLinesOn             = False         ; turn off contour lines
+  resmap@cnLevelSelectionMode  = "ManualLevels"; Manually set contour levels.
+  resmap@cnMinLevelValF        =  5.
+;  resmap@cnMaxLevelValF               = 3200.
+;  resmap@cnLevelSpacingF              =  200.
+
+;=================================================;
+; PLOT ZS 
+;=================================================;
+  plot_zs = gsn_csm_contour_map(wks,ZS,resmap)
+
+end
+
+
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_Cloud.ncl b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_Cloud.ncl
new file mode 100644 (file)
index 0000000..ac7a1e4
--- /dev/null
@@ -0,0 +1,217 @@
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+
+; LCARTESIAN = False
+
+begin
+;=================================================;
+; Open file
+;=================================================;
+; MESONH input file
+  mnh_file="REUNI.1.00A20.004.nc4"
+  a = addfile(mnh_file, "r")
+
+;==================================================;
+; Open the workstation and choose colormap
+; For paper quality plot do not use type ncgm or eps
+; Use type= ps or pdf 
+;==================================================;
+
+  type = "x11"          ; could be png, ps , pdf ...
+  wks = gsn_open_wks(type,"plt_Cloud")
+
+; Choose the colormap you want 
+; See http://www.ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml
+; for all the available color tables
+  gsn_define_colormap(wks,"WhBlGrYeRe")   ; overwrite the .hluresfile color map
+
+;=================================================;
+; Get informations on variable sizes
+;=================================================;
+
+  mdims = getfilevardimsizes(a,"UT") ; get some dimension sizes for the file
+  nd = dimsizes(mdims)
+  imax=mdims(nd-1)-2
+  jmax=mdims(nd-2)-2
+  kmax=mdims(nd-3)-2
+; -2 : to remove non-physical values
+
+;=================================================;
+; Read the variables we will need
+;=================================================;
+
+; Read latitude and longitude fields
+; -----------------------------------------
+  lat2d = a->LAT(1:jmax,1:imax)
+  lat2d@units="degrees_north"
+  lon2d = a->LON(1:jmax,1:imax)
+  lon2d@units="degrees_east"
+
+; Read water wapor mixing ratio if it exists
+; -----------------------------------------
+    if(isfilevar(a,"RVT"))
+      qv = a->RVT(1:kmax,1:jmax,1:imax)
+      qv = qv*1000.
+      qv@units = "g/kg"   
+      qv@description="Water vapor mixing ratio"
+      qv@lat2d=lat2d
+      qv@lon2d=lon2d
+    end if
+
+; Read cloud mixing ratio if it exists
+; -----------------------------------------
+    if(isfilevar(a,"RCT"))
+      qc = a->RCT(1:kmax,1:jmax,1:imax)
+      qc = qc*1000.
+      qc@units = "g/kg"   
+      qc@description="Cloud mixing ratio"
+      qc@lat2d=lat2d
+      qc@lon2d=lon2d
+    end if
+
+; Read rain mixing ratio if it exists
+; -----------------------------------------
+    if(isfilevar(a,"RRT"))
+      qr = a->RRT(1:kmax,1:jmax,1:imax)
+      qr = qr*1000.
+      qr@units = "g/kg"   
+      qr@description="Rain mixing ratio"
+      qr@lat2d=lat2d
+      qr@lon2d=lon2d
+    end if
+
+; Read ice mixing ratio if it exists
+; -----------------------------------------
+    if(isfilevar(a,"RIT"))
+      qi = a->RIT(1:kmax,1:jmax,1:imax)
+      qi = qi*1000.
+      qi@units = "g/kg"   
+      qi@description="Ice mixing ratio"
+      qi@lat2d=lat2d
+      qi@lon2d=lon2d
+    end if
+
+; Read projection parameters
+; --------------------
+    RPK  = a->RPK
+    BETA = a->BETA
+    LON0 = a->LON0
+
+; Set resource for the map projection
+;-----------------------------------------------;
+  resmap=True
+
+; Set map projection ressources using projection parameters
+;-----------------------------------------------;
+  if (RPK.gt.0)
+; ---------------------------
+  ;   Lambert  projection from north pole
+; ---------------------------
+   resmap@mpProjection          = "LambertConformal"     ; projection
+   pole                         = 1    ; projection for north hemisphere
+   resmap@mpLambertParallel1F   = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere
+   resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file
+   resmap@mpLambertMeridianF    = LON0      ; ncl adds from grib file
+  end if
+
+  if (RPK.lt.0)
+; ---------------------------
+  ;   Lambert projection from south pole
+; ---------------------------
+   resmap@mpProjection          = "LambertConformal"     ; projection
+   pole                         = -1                     ; projection for south hemisphere
+   resmap@mpLambertParallel1F   = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere
+   resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file
+   resmap@mpLambertMeridianF    =  LON0      ; ncl adds from grib file
+  end if
+
+  if (RPK.eq.1)
+; ---------------------------
+  ;   Stereographic projection north
+; ---------------------------
+    resmap@mpProjection = "Stereographic"
+    resmap@mpCenterLonF           = LON0
+    resmap@mpCenterRotF           = BETA
+    resmap@mpCenterLatF           = 90
+  end if
+
+  if (RPK.eq.-1)
+; ---------------------------
+  ;   Stereographic projection south
+; ---------------------------
+    resmap@mpProjection = "Stereographic"
+    resmap@mpCenterLonF           = LON0
+    resmap@mpCenterRotF           = BETA
+    resmap@mpCenterLatF           = -90
+  end if
+
+  if (RPK.eq.0) then
+; ---------------------------
+  ;   Mercator projection
+; ---------------------------
+    resmap@mpProjection = "Mercator"
+  end if
+
+ print("Map projection="+resmap@mpProjection)
+
+;=================================================;
+; Set some other basic resources
+;=================================================;
+
+  resmap@mpDataBaseVersion      = "HighRes"  ; choose highres map data version (must be donwloaded)
+  resmap@mpFillOn               = False
+
+; Defining the corners for projection
+;====================================
+  resmap@mpLimitMode            = "Corners"
+  resmap@mpLeftCornerLatF       = lat2d(0,0)
+  resmap@mpLeftCornerLonF       = lon2d(0,0)
+  resmap@mpRightCornerLatF     = lat2d(jmax-1,imax-1)
+  resmap@mpRightCornerLonF     = lon2d(jmax-1,imax-1)
+
+
+;=========================================================;
+;  Loop over levels by step of 5
+;=========================================================;
+    do level = 0,19,5      
+
+      display_level = level + 1
+
+; Contour resources
+;--------------------
+      opts           =  resmap
+      opts@cnFillOn  = True
+      opts@cnLinesOn = False
+
+      opts@tiMainString = "Model Level "+display_level
+
+
+; Plot wapor mixing ratio if it exists
+; -----------------------------------------
+      if (isvar("qv"))
+        plot=gsn_csm_contour_map(wks,qv(level,:,:),opts)
+      end if
+
+; Plot cloud mixing ratio if it exists
+; -----------------------------------------
+      if (isvar("qc"))
+        plot=gsn_csm_contour_map(wks,qc(level,:,:),opts)
+      end if
+
+; Plot rain mixing ratio if it exists
+; -----------------------------------------
+      if (isvar("qr"))
+        plot=gsn_csm_contour_map(wks,qr(level,:,:),opts)
+      end if
+
+; Plot ice mixing ratio if it exists
+; -----------------------------------------
+      if (isvar("qi"))
+        plot=gsn_csm_contour_map(wks,qi(level,:,:),opts)
+      end if
+
+
+    end do      ; END OF LEVEL LOOP
+;=================================================;
+
+end
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_CrossSection.ncl b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_CrossSection.ncl
new file mode 100644 (file)
index 0000000..938c64d
--- /dev/null
@@ -0,0 +1,314 @@
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "MESONHtools.ncl"
+
+begin
+
+;=================================================;
+; Open file 
+;=================================================;
+  in_file="REUNI.1.00A20.004.nc4"
+  a = addfile(in_file, "r")
+
+;=================================================;
+; We generate plots, but what kind do we prefer?
+; For paper quality plot do not use ncgm or eps
+;=================================================;
+  type = "x11"
+;  type = "png"
+; type = "pdf"
+; type = "ps"
+
+;=================================================;
+; Open the workstation and choose colormap
+;=================================================;
+  wks = gsn_open_wks(type,"plt_CrossSection2")
+  gsn_define_colormap(wks,"topo_15lev")
+
+;=================================================;
+; Get informations on variable sizes
+;=================================================;
+  mdims = getfilevardimsizes(a,"UT") ; get some dimension sizes for the file
+  nd = dimsizes(mdims)
+  imax=mdims(nd-1)-2
+  jmax=mdims(nd-2)-2
+  kmax=mdims(nd-3)-2
+
+;=================================================;
+; Read  some variables from mnh_file
+;=================================================;
+    ZS =a->ZS(1:jmax,1:imax)              
+
+    th = a->THT(1:kmax,1:jmax,1:imax)
+    th@long_name="Potential Temperature"
+    th@description="TH"
+
+    um = a->UT(1:kmax,1:jmax,1:imax+1)
+    um@description="U"
+
+    zh= a->ZHAT(0:kmax)
+
+    lat2d = a->LAT(1:jmax,1:imax)
+    lat2d@units="degrees_north"
+    lon2d = a->LON(1:jmax,1:imax)
+    lon2d@units="degrees_east"
+
+    ZS@lat2d=lat2d
+    ZS@lon2d=lon2d
+
+    th@lat2d=lat2d
+    th@lon2d=lon2d
+
+    um@lat2d=lat2d
+    um@lon2d=lon2d
+
+;=================================================;
+; Take care of the Arakawa C grid
+;=================================================;
+; Unstagger zh (from grid 4 to 1)
+;----------------------------------;
+    nzh=new(kmax,double)
+    do k=0,kmax-1
+     nzh(k)=(zh(k)+zh(k+1))/2.
+    end do
+
+; Unstagger U
+;----------------------------------;
+  um!0="Z"
+  um!1="Y"
+  um!2="X"
+  u_unst = wrf_user_unstagger(um,"X")
+
+;=================================================;
+; Needed for the wrf_intrp3d function
+;=================================================;
+; Create Z3D 
+;----------------------------------;
+    z=new(dimsizes(th),double)
+    zcoef=1.-ZS/nzh(kmax-1)
+
+    do i=0,imax-1
+    do j=0,jmax-1
+       z(:,j,i) = nzh*zcoef(j,i)+ZS(j,i)
+    end do
+    end do
+
+; get height info for labels
+;----------------------------------;
+      zmin = 0.
+      zmax = 6.         ; We are only interested in the first 6km
+      nz   = floattoint(zmax + 1)
+
+
+;----------------------------------------------------------;
+
+; Resources for the map projection
+;-----------------------------------------------;
+  resmap=True
+
+;=================================================;
+; Define Other Resources for map ploting
+;=================================================;
+
+  resmap@gsnFrame               = False
+  resmap@gsnMaximize            = True
+
+; Map resources
+;-----------------------------------------------;
+  resmap@mpFillOn               = False
+  resmap@mpGridAndLimbOn        = True            ; turn on lat/lon lines
+  resmap@mpGridLatSpacingF      = 10.             ; spacing for lat lines
+  resmap@mpGridLonSpacingF      = 10.             ; spacing for lon lines
+  resmap@mpDataBaseVersion      = "HighRes"     ; choose highres map data version (must be donwloaded)
+
+; Defining the corners for projection
+;====================================
+  resmap@mpLimitMode           = "Corners"
+  resmap@mpLeftCornerLatF      = lat2d(0,0)
+  resmap@mpLeftCornerLonF      = lon2d(0,0)
+  resmap@mpRightCornerLatF     = lat2d(jmax-1,imax-1)
+  resmap@mpRightCornerLonF     = lon2d(jmax-1,imax-1)
+
+
+; Contour definitions
+;-----------------------------------------------;
+  resmap@cnFillOn               = True
+  resmap@cnLinesOn              = False
+  resmap@cnLevelSelectionMode   = "ExplicitLevels"
+  cnLevels                      = (/1,250,500,750,1000,1250,1500,1750,2000,2250,2500,2750,3000/)
+  resmap@cnLevels               = cnLevels
+
+  resmap@lbOrientation          = "Vertical"
+
+  
+  FirstPlot = True
+
+;=================================================;
+; Planes definition and interpolation
+;=================================================;
+    do ip = 1, 3       ; we are doing 3 plots
+                       ; all with the pivot point (plane) in the center of the domain
+                       ; at angles 0, 45 and 90
+ ;  
+ ;                   |
+ ;       angle=0 is  |
+ ;                   |
+ ; 
+
+; Define the plane where the cross-section will be done
+; ------------------------------------------------------
+        plane = new(2,float)
+        plane = (/ (mdims(nd-1)-2)/2, (mdims(nd-2)-2)/2 /)    ; pivot point is center of domain (x,y)
+
+; Resource for the interpolation
+; --------------------------------
+        opts = False
+
+        if(ip .eq. 1) then
+          angle = 90.
+          X_plane = wrf_user_intrp2d(lon2d,plane,angle,opts)
+          X_desc = "longitude"
+        end if
+        if(ip .eq. 2) then
+          angle = 0.
+          X_plane = wrf_user_intrp2d(lat2d,plane,angle,opts)
+          X_desc = "latitude"
+        end if
+        if(ip .eq. 3) then
+          angle = -45.
+          X_plane = wrf_user_intrp2d(lon2d,plane,angle,opts)
+          X_desc = "longitude"
+          
+        end if
+
+; Interpolate fields onto the defined plane
+; --------------------------------------------
+        um_plane = wrf_user_intrp3d(u_unst,z,"v",plane,angle,opts)
+        th_plane = wrf_user_intrp3d(th,z,"v",plane,angle,opts)
+
+; Find the index where 6km is - only need to do this once
+; --------------------------------------------
+        if ( FirstPlot ) then
+          zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)
+          b = ind(zz(:,0) .gt. zmax*1000. )
+          zmax_pos = b(0) - 1
+          if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then
+            zspan = b(0) - 1
+          else
+            zspan = b(0)
+          end if
+          delete(zz)
+          delete(b)
+          FirstPlot = False
+        end if
+
+; Define X-axis labels
+;----------------------
+      dimsX = dimsizes(X_plane)
+      xmin  = X_plane(0)
+      xmax  = X_plane(dimsX(0)-1)
+      xspan = dimsX(0)-1
+      nx    = 3
+
+; Options for XY Plots
+;-----------------------------------------
+        opts_xy                         = True
+        opts_xy@tiXAxisString           = X_desc
+        opts_xy@tiYAxisString           = "Height (km)"
+        opts_xy@tiXAxisFontHeightF      = 0.020
+        opts_xy@tiYAxisFontHeightF      = 0.020
+
+; Resources to plot the topography (missing values in fields)
+;-----------------------------------------
+        opts_xy@cnMissingValPerimOn     = True
+        opts_xy@cnMissingValFillColor   = "red"
+        opts_xy@cnMissingValFillPattern = 8
+    
+; Tickmarks/labels resources
+;------------------------------
+        opts_xy@lbOrientation           = "Vertical"
+
+        opts_xy@tmXTOn                  = False
+        opts_xy@tmYROn                  = False
+        opts_xy@tmXBMode                = "Explicit"
+        opts_xy@tmXBValues              = fspan(0,xspan,nx)              ; Create tick marks
+        opts_xy@tmXBLabels              = sprintf("%.2f",fspan(xmin,xmax,nx))  ; Create labels
+        opts_xy@tmXBLabelFontHeightF    = 0.015
+        opts_xy@tmYLMode                = "Explicit"
+        opts_xy@tmYLValues              = fspan(0,zspan,nz)                    ; Create tick marks
+        opts_xy@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz))  ; Create labels
+        opts_xy@tmXBMajorLengthF        = 0.02
+        opts_xy@tmYLMajorLengthF        = 0.02
+        opts_xy@tmYLLabelFontHeightF    = 0.015
+
+
+; Plotting options for U
+;--------------------------
+        opts_um                         = opts_xy
+        opts_um@cnFillOn                = True
+
+; Plotting options for Temperature
+;--------------------------
+        opts_th                           = opts_xy
+        opts_th@cnInfoLabelZone           = 1
+        opts_th@cnInfoLabelSide           = "Top"
+        opts_th@cnInfoLabelPerimOn        = True
+        opts_th@cnInfoLabelOrthogonalPosF = -0.00005
+
+        opts_th@gsnLeftStringParallelPosF = 0.1
+        opts_th@gsnLeftStringOrthogonalPosF = 0.01
+
+;======================
+; MAKE PLOTS         
+;======================
+
+; Horizontal plot with topography
+;-----------------------------------------------------------------------
+        plot=gsn_csm_contour_map(wks,ZS,resmap)
+
+; Draw line that shows the position of the cross section
+;-----------------------------------------------------------------------
+        lat_plane = wrf_user_intrp2d(lat2d,plane,angle,opts)
+        lon_plane = wrf_user_intrp2d(lon2d,plane,angle,opts)
+
+        lnres = True
+        lnres@gsLineThicknessF = 3.0
+        lnres@gsLineColor = "Red"
+        do ii = 0,dimsX(0)-2
+            gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres)
+        end do
+        frame(wks)
+        delete(lon_plane)
+        delete(lat_plane)
+
+; Vertical cross-section
+;--------------------------
+        opts_th@gsnFrame = False
+        opts_th@gsnDraw  = False
+
+        
+
+        opts_um@gsnFrame = False
+        opts_um@gsnDraw  = False
+
+        contour_th = gsn_csm_contour(wks,th_plane(0:zmax_pos,:),opts_th)
+        contour_um = gsn_csm_contour(wks,um_plane(0:zmax_pos,:),opts_um)
+
+        overlay(contour_um,contour_th)    ; plot x-section
+        draw(contour_um)
+        frame(wks)
+
+
+; Delete options and fields, so we don't have carry over
+;-------------------------------------------------------
+        delete(opts_xy)
+        delete(opts_th)
+        delete(opts_um)
+        delete(th_plane)
+        delete(um_plane)
+        delete(X_plane)
+
+end do  ; make next cross section
+
+
+end
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_ModelLevels.ncl b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_ModelLevels.ncl
new file mode 100644 (file)
index 0000000..49640ea
--- /dev/null
@@ -0,0 +1,320 @@
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "MESONHtools.ncl"
+
+;=============================================================;
+;  Example of script to plot field on model levels from MESONH 
+;  netcdf file generated with lfi2cdf
+;=============================================================;
+
+begin
+
+;=================================================;
+; Open file 
+;=================================================;
+; The MESONH input file.  
+; ----------------------
+  mnh_file="REUNI.1.00A20.004.nc4"
+  a = addfile(mnh_file,"r")
+
+;==================================================;
+; Open the workstation and choose colormap
+; For paper quality plot do not use type ncgm or eps
+; Use ps or pdf or x11 (for debugging)
+;==================================================;
+  type = "x11"
+  wks = gsn_open_wks(type,"plt_ModelLevels")
+  gsn_define_colormap(wks,"rainbow+gray")
+
+;=================================================;
+; Get informations on variable sizes
+; dims are dims-2 to remove non-physical values
+;=================================================;
+  mdims = getfilevardimsizes(a,"THT") ; get dimension sizes
+  nd = dimsizes(mdims)
+  imax=mdims(nd-1)-2
+  jmax=mdims(nd-2)-2
+  kmax=mdims(nd-3)-2
+
+;=================================================;
+; Read the variables we need
+;=================================================;
+  lat2d = a->LAT(1:jmax,1:imax)
+  lat2d@units="degrees_north"
+  lon2d = a->LON(1:jmax,1:imax)
+  lon2d@units="degrees_east"
+
+  th  = a->THT(1:kmax,1:jmax,1:imax)   ; theta
+  th@long_name="Potential Temperature"
+  th@units = "K"
+  th@lat2d=lat2d
+  th@lon2d=lon2d
+
+  qv  = a->RVT(1:kmax,1:jmax,1:imax)   ; Qv
+  qv = qv*1000.
+  qv@long_name="Water vapor mixing ratio"
+  qv@units = "g/kg"
+  qv@lat2d=lat2d
+  qv@lon2d=lon2d
+
+  ut   = a->UT(1:kmax,1:jmax,1:imax+1)        ; u
+  ut@long_name="U"
+  ut@units="m/s"
+  ut@lat2d=lat2d
+  ut@lon2d=lon2d
+
+  vt   = a->VT(1:kmax,1:jmax+1,1:imax)        ; v
+  vt@long_name="V"
+  vt@units="m/s"
+  vt@lat2d=lat2d
+  vt@lon2d=lon2d
+
+ ; Unstagger U
+; --------------------
+    ut!0="Z"
+    ut!1="Y"
+    ut!2="X"
+    u = wrf_user_unstagger(ut,"X")
+
+ ; Unstagger V
+; --------------------
+    vt!0="Z"
+    vt!1="Y"
+    vt!2="X"
+    v = wrf_user_unstagger(vt,"Y")
+    v@description="V"
+
+; Calculate wind speed
+; --------------------
+    spd = (u*u + v*v)^(0.5)             ; speed in m/sec
+    spd@long_name = "Wind Speed"
+    spd@units = "m/s"
+    spd@lat2d = lat2d
+    spd@lon2d = lon2d
+
+; Read projection parameters
+; --------------------
+    RPK  = a->RPK
+    BETA = a->BETA
+    LON0 = a->LON0
+
+
+; Set resource for the map projection
+;-----------------------------------------------;
+  resmap=True
+
+; Set map projection ressources using projection parameters
+;-----------------------------------------------;
+  if (RPK.gt.0)
+; ---------------------------
+  ;   Lambert  projection from north pole
+; ---------------------------
+   resmap@mpProjection          = "LambertConformal"     ; projection
+   pole                         = 1    ; projection for north hemisphere
+   resmap@mpLambertParallel1F   = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere
+   resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file
+   resmap@mpLambertMeridianF    = LON0      ; ncl adds from grib file
+  end if
+
+  if (RPK.lt.0)
+; ---------------------------
+  ;   Lambert projection from south pole
+; ---------------------------
+   resmap@mpProjection          = "LambertConformal"     ; projection
+   pole                         = -1                     ; projection for south hemisphere
+   resmap@mpLambertParallel1F   = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere
+   resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file
+   resmap@mpLambertMeridianF    =  LON0      ; ncl adds from grib file
+  end if
+
+  if (RPK.eq.1)
+; ---------------------------
+  ;   Stereographic projection north
+; ---------------------------
+    resmap@mpProjection = "Stereographic"
+    resmap@mpCenterLonF           = LON0
+    resmap@mpCenterRotF           = BETA
+    resmap@mpCenterLatF           = 90
+  end if
+
+  if (RPK.eq.-1)
+; ---------------------------
+  ;   Stereographic projection south
+; ---------------------------
+    resmap@mpProjection = "Stereographic"
+    resmap@mpCenterLonF           = LON0
+    resmap@mpCenterRotF           = BETA
+    resmap@mpCenterLatF           = -90
+  end if
+
+  if (RPK.eq.0) then
+; ---------------------------
+  ;   Mercator projection
+; ---------------------------
+    resmap@mpProjection = "Mercator"
+  end if
+
+ print("Map projection="+resmap@mpProjection)
+
+; Defining the corners for projection
+; --------------------------------
+  resmap@mpLimitMode            = "Corners"
+  resmap@mpLeftCornerLatF       = lat2d(0,0)
+  resmap@mpLeftCornerLonF       = lon2d(0,0)
+  resmap@mpRightCornerLatF     = lat2d(jmax-1,imax-1)
+  resmap@mpRightCornerLonF     = lon2d(jmax-1,imax-1)
+
+;=================================================;
+; Set some other basic resources
+;=================================================;
+  resmap@mpDataBaseVersion      = "HighRes"  ; highres map data version 
+  resmap@mpFillOn               = False
+  resmap@mpGeophysicalLineColor = "Black"
+  resmap@mpNationalLineColor    = "Black"
+
+  resmap@pmTickMarkDisplayMode = "Always"         ; turn on tickmarks
+
+;=========================================================;
+;  Loop over levels by step of 20
+;=========================================================;
+    do level =0,kmax-1,20      ; LOOP OVER LEVELS
+
+       display_level          = level + 1
+       resmap@tiMainString    = "Model Level  " + display_level
+
+;=======================
+; Resources for Theta
+;=======================
+       opts                                = resmap
+
+       opts@gsnFrame              = False
+       opts@gsnDraw               = False
+       opts@gsnContourLineThicknessesScale = 2.0
+
+       opts@cnLineColor                    = "Red"
+       opts@cnInfoLabelOrthogonalPosF      = -0.05  
+       opts@cnInfoLabelParallelPosF        = 0.97  
+
+       contour = gsn_csm_contour_map(wks,th(level,:,:),opts)
+       draw(contour)
+       frame(wks)
+       delete(opts)
+  
+;======================
+; Resources for Qv
+;======================
+       opts = resmap
+
+; General resources
+;-------------------
+       opts@gsnFrame              = False
+       opts@gsnDraw               = False
+       opts@gsnSpreadColors       = True
+       opts@gsnSpreadColorStart   = 30
+       opts@gsnSpreadColorEnd     = -2
+
+; Contour resources
+;-------------------
+       opts@cnLinesOn         = False   ; no line
+       opts@cnFillOn          = True    ; filled
+
+       contour = gsn_csm_contour_map(wks,qv(level,:,:),opts)
+       draw(contour)
+       frame(wks)
+       delete(opts)
+   
+;===============================================
+; Overlay Wind Vectors and Speed over basic map
+;===============================================
+
+;=======================
+; 1.Create Base map plot
+;=======================
+  resmap@gsnDraw                = False         ; don't draw yet
+  resmap@gsnFrame               = False         ; don't advance frame yet
+
+  map=gsn_map(wks,"Stereographic",resmap)
+
+;=======================
+; 2.Create Wind speed plot
+;=======================
+  opts_ws = True
+
+; General resources
+;-------------------
+  opts_ws@gsnFrame              = False
+  opts_ws@gsnDraw               = False
+  opts_ws@gsnSpreadColors       = True
+  opts_ws@gsnSpreadColorStart   = 30
+  opts_ws@gsnSpreadColorEnd     = -2
+
+; Tickmark resources
+;---------------------
+; No tickmark and label on axis for this plot
+; They are already drawn by the basic map plot
+; -------------------------------------------
+;  opts_ws@tmXBOn            = False
+;  opts_ws@tmYROn            = False
+;  opts_ws@tmXBLabelsOn      = False
+;  opts_ws@tmYRLabelsOn      = False
+
+  opts_ws@lbOrientation     = "Vertical"
+
+
+; Transformation resources
+;---------------------------
+; Field is already on the right projection
+;------------------------------------------
+;  opts_ws@tfDoNDCOverlay = True
+
+; Contour resources
+; -------------------
+  opts_ws@cnFillOn     = True
+  opts_ws@cnLinesOn    = False
+
+; PLOT
+;-----
+  contour = gsn_csm_contour(wks,spd(level,:,:),opts_ws)
+
+;===========================
+; 3.Create Wind vector plot
+;===========================
+  opts_vec = True
+
+  opts_vec@gsnFrame   = False
+  opts_vec@gsnDraw    = False
+  opts_vec@gsnLeftString = "Wind Vector"
+  opts_vec@gsnLeftStringParallelPosF = 0.3 
+  opts_vec@gsnLeftStringOrthogonalPosF = 0.018 
+
+; No Tickmarks/labels
+;-------------
+;  opts_vec@tmXBOn            = False
+;  opts_vec@tmYROn            = False
+;  opts_vec@tmXBLabelsOn      = False
+;  opts_vec@tmYRLabelsOn      = False
+
+;  opts_vec@tfDoNDCOverlay = True
+
+; Vector resources
+;-------------------
+  opts_vec@vcRefLengthF      = 0.1
+  opts_vec@vcRefMagnitudeF   = 20
+  opts_vec@vcMinDistanceF    = 0.05
+
+
+; PLOT
+;-----
+   vector =  gsn_csm_vector(wks,u(level,:,:),v(level,:,:),opts_vec)
+
+;===================
+; 4.Overlay and draw
+;===================
+   overlay(map,contour)
+   overlay(map,vector)
+
+   draw(map)
+   frame(wks)
+
+ end do      ; END OF LEVEL LOOP
+   
+end
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_PressureLevel.ncl b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/mnh_PressureLevel.ncl
new file mode 100644 (file)
index 0000000..a992ef2
--- /dev/null
@@ -0,0 +1,246 @@
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "MESONHtools.ncl"
+
+;=================================================;
+; Example of script to interpolate from model to pressure level 
+; and to plot some fields from MESONH netcdf4 files
+;=================================================;
+
+begin
+
+;=======================================================;
+; Pressure levels that we want the data interpolated to
+;=======================================================;
+  pressure_levels = (/ 850., 700., 500., 300./)  
+  nlevels         = dimsizes(pressure_levels)  ; number of pressure levels
+
+;=================================================;
+; Open file 
+; MESONH netcdf4 input file 
+;=================================================;
+
+  mnh_file="REUNI.1.00A20.004.nc4"
+  a = addfile(mnh_file, "r")
+
+;=================================================;
+; Open the workstation and choose colormap
+; For paper quality plot do not use type ncgm or eps
+; Use ps or pdf and x11 (for debugging)
+;==================================================;
+
+  type = "x11"
+  wks = gsn_open_wks(type,"plt_PressureLevel")
+  gsn_define_colormap(wks,"spread_15lev")
+
+;=================================================;
+; Get informations on variable sizes
+; dims are dims-2 to remove non-physical values
+;=================================================;
+  print("Reading dims")
+  mdims = getfilevardimsizes(a,"UT") ; get dimension sizes
+  nd = dimsizes(mdims)
+  imax=mdims(nd-1)-2
+  jmax=mdims(nd-2)-2
+  kmax=mdims(nd-3)-2
+
+
+;=================================================;
+; Read the variables we need
+;=================================================;
+
+  print("Reading THT")
+  tk = a->THT(1:kmax,1:jmax,1:imax)         ; potential temperature
+  tk@description="Potential Temperature"
+
+  print("Reading wind")
+  u = a->UT(1:kmax,1:jmax,1:imax+1)         ; u
+  v = a->VT(1:kmax,1:jmax+1,1:imax)         ; v
+
+; Unstagger u,v
+  u!0="Z"
+  u!1="Y"
+  u!2="X"
+  u_unst = wrf_user_unstagger(u,"X")
+
+  v!0="Z"
+  v!1="Y"
+  v!2="X"
+  v_unst = wrf_user_unstagger(v,"Y")
+
+  print("Reading pressure")
+  p = a->PABST(1:kmax,1:jmax,1:imax)        ; pressure
+  p@description="Pressure"
+
+  print("Reading water vapor mixing ratio")
+  rv = a->RVT(1:kmax,1:jmax,1:imax)    ; water vapor mixing ratio
+  rv = rv*1000.    ; water vapor mixing ratio
+  rv@description="Water vapor mixing ratio"
+
+  print("Reading orography")
+  zs = a->ZS(1:jmax,1:imax)                ; terrain
+  zh = a->ZHAT(1:kmax+1)                   ; heigth without terrain
+
+ print("End of variables reading")
+
+;==========================================
+; Create 3D Z to be used for interpolation
+;==========================================
+
+  z=new(dimsizes(tk),double)
+  z@description="Height"
+
+; Unstag zh
+;----------
+   nzh=new(kmax,double)
+   do k=0,kmax-1
+   nzh(k)=(zh(k)+zh(k+1))/2.
+   end do
+   zcoef=1.-zs/nzh(kmax-1)
+
+    do i=0,imax-1
+    do j=0,jmax-1
+       z(:,j,i) = nzh*zcoef(j,i)+zs(j,i)
+    end do
+    end do
+
+;=================================================;
+; Horizontal interpolation on pressure levels
+;=================================================;
+
+
+    do level = 0,nlevels-1                 ; LOOP OVER LEVELS
+
+      pressure = pressure_levels(level)
+
+      tk_plane = wrf_user_intrp3d(tk,p,"h",pressure,0.,False)
+      z_plane  = wrf_user_intrp3d( z,p,"h",pressure,0.,False)
+      u_plane  = wrf_user_intrp3d( u_unst,p,"h",pressure,0.,False)
+      v_plane  = wrf_user_intrp3d( v_unst,p,"h",pressure,0.,False)
+      rv_plane  = wrf_user_intrp3d( rv,p,"h",pressure,0.,False)
+
+      spd_plane = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec
+      spd_plane@description = "Wind Speed"
+      spd_plane@units = "m/s"
+
+      u_plane@units = "m/s"
+      v_plane@units = "m/s"
+
+;==================================
+; Create plots usinf wrf functions
+;==================================
+      res=True
+      res@MainTitle                   = "MESONH"
+      res@Footer = False
+
+; Plotting options for Tk                
+;--------------------------
+      opts =  res                         
+      opts@cnLineColor = "Red"
+      opts@cnInfoLabelOrthogonalPosF = 0.07 ; offset second label information
+      opts@gsnContourLineThicknessesScale = 2.0
+
+      contour_tk = wrf_contour(a,wks,tk_plane,opts)
+      delete(opts)
+
+; Plotting options for  RV          
+;--------------------------
+  opts                          = res                          
+  opts@cnFillOn                 = True  
+  opts@pmLabelBarOrthogonalPosF = -0.1
+  opts@cnMissingValFillColor    = "Red"
+  opts@cnMissingValFillPattern  = 8
+
+  contour_rv = wrf_contour(a,wks,rv_plane,opts)
+  delete(opts)
+
+; Plotting options for Wind Speed                
+;--------------------------
+  opts = res                          
+  opts@cnLineColor = "MediumSeaGreen"
+  opts@cnInfoLabelOrthogonalPosF = 0.07  ; offset second label information
+  opts@gsnContourLineThicknessesScale = 3.0
+
+  contour_spd = wrf_contour(a,wks,spd_plane,opts)
+  delete(opts)
+
+
+; Plotting options for Wind Vectors                 
+;-----------------------------------
+  opts            = res          
+  opts@FieldTitle = "Wind"   ; overwrite Field Title
+  opts@NumVectors = 47       ; wind barb density
+
+  vector = wrf_vector(a,wks,u_plane,v_plane,opts)
+  delete(opts)
+
+; Plotting options for Geopotential Heigh
+;-----------------------------------------
+  opts_z                                = res                          
+  opts_z@cnLineColor                    = "Blue"
+  opts_z@gsnContourLineThicknessesScale = 2.0
+
+;=============
+; MAKE PLOTS                                       
+;=============
+
+  pltres = True
+
+; Resources for map
+;-------------------
+  resmap                             = True
+  resmap@mpGeophysicalLineColor      = "Black"
+  resmap@mpNationalLineColor         = "Black"
+  resmap@mpGridLineColor             = "Black"
+  resmap@mpLimbLineColor             = "Black"
+  resmap@mpPerimLineColor            = "Black"
+  resmap@mpGeophysicalLineThicknessF = 2.0
+  resmap@mpGridLineThicknessF        = 2.0
+  resmap@mpLimbLineThicknessF        = 2.0
+  resmap@mpNationalLineThicknessF    = 2.0
+  resmap@mpGridAndLimbOn             = False            ; turn off lat/lon lines
+  resmap@mpDataBaseVersion           = "HighRes"     ; choose highres map data
+
+
+; Overlay rv , tk , height and wind barbs
+;------------------------------------------
+  if ( pressure .eq. 850 ) then   ; plot temp, rv, height, wind barbs
+
+   contour_height = wrf_contour(a,wks,z_plane,opts_z)
+   plot = mnh_map_overlays(a,wks,(/contour_rv,contour_tk,contour_height, \
+                                    vector/),pltres,resmap)
+  end if
+
+; Overlay  tk , height and wind barbs
+;------------------------------------------
+  if ( pressure .eq. 700 ) then   ; plot temp, height, wind barbs
+
+   contour_height = wrf_contour(a,wks, z_plane,opts_z)
+   plot = mnh_map_overlays(a,wks,(/contour_tk,contour_height, \
+                                    vector/),pltres,resmap)
+  end if
+
+; Overlay  tk , height and wind barbs
+;------------------------------------------
+  if ( pressure .eq. 500 ) then   ; plot temp, height, wind barbs
+
+   contour_height = wrf_contour(a,wks, z_plane,opts_z)
+   plot = mnh_map_overlays(a,wks,(/contour_tk,contour_height, \
+                                    vector/),pltres,resmap)
+  end if
+
+; Overlay  wind speed , height and wind barbs
+;------------------------------------------
+  if ( pressure .eq. 300 ) then   ; plot windspeed, height, wind barbs
+
+   contour_height = wrf_contour(a,wks, z_plane,opts_z)
+   plot = mnh_map_overlays(a,wks,(/contour_spd,contour_height, \
+                                    vector/),pltres,resmap)
+  end if
+
+  delete(opts_z)
+
+end do      ; END OF LEVEL LOOP
+
+
+end
diff --git a/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/run_ncl b/MY_RUN/KTEST/004_Reunion/005_ncl_nc4/run_ncl
new file mode 100755 (executable)
index 0000000..7ed60c1
--- /dev/null
@@ -0,0 +1,26 @@
+set -x
+#### Get highres coastline data for ncl
+# If you don't have internet connexion, comment the next line
+# And change in the ncl script HighRes by MediumRes
+./get_ncl_highres_files
+#### nc4 file to plot
+ln -sf ../003_mesonh/REUNI.1.00A20.004.nc4 .
+#### Initilizing variables for ncl
+path_ncl=`which ncl`
+if [ "$path_ncl" == "" ]
+then
+echo "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
+echo "!   ERROR: NCL is not installed or the path to   !"
+echo "!     ncl binary is not set correctly            !"
+echo "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
+exit
+fi
+export NCARG_ROOT=${path_ncl%%/bin/ncl}
+export PATH=$NCARG_ROOT/bin:$PATH
+export NCARG_RANGS=.
+#### Running  ncl scripts 
+ncl mnh_BasicMap.ncl
+ncl mnh_ModelLevels.ncl
+ncl mnh_Cloud.ncl
+ncl mnh_CrossSection.ncl
+ncl mnh_PressureLevel.ncl
diff --git a/MY_RUN/KTEST/004_Reunion/006_ncl/clean_ncl b/MY_RUN/KTEST/004_Reunion/006_ncl/clean_ncl
new file mode 100755 (executable)
index 0000000..7b1634b
--- /dev/null
@@ -0,0 +1,2 @@
+
+rm -f  *.nc4 *.ps 
diff --git a/MY_RUN/KTEST/004_Reunion/006_ncl/plot_Reunion.ncl b/MY_RUN/KTEST/004_Reunion/006_ncl/plot_Reunion.ncl
new file mode 100644 (file)
index 0000000..9898f34
--- /dev/null
@@ -0,0 +1,320 @@
+;================================================;
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
+; ================================================;
+;================================================;
+  begin
+;=================================================;
+; open file and read in data
+;=================================================;
+  fichier1 = addfile("REUNI.1.00A20.004.nc4", "r")
+;==================================================;
+; Open the workstation
+;==================================================;
+  type = "ps"
+  wks = gsn_open_wks(type,"visu_Reunion")
+;=================================================;
+; Get informations on variable sizes
+; dims are dims-2 to remove non-physical values
+;=================================================;
+  mdims = getfilevardimsizes(fichier1,"UT") ; get dimension sizes
+  nd = dimsizes(mdims)
+  imax=mdims(nd-1)-2
+  jmax=mdims(nd-2)-2
+  kmax=mdims(nd-3)-2
+;=================================================;
+; Set map projection
+;=================================================;
+  lat2d=fichier1->LAT(1:jmax,1:imax)
+  lon2d=fichier1->LON(1:jmax,1:imax)
+
+; Resources for the map projection
+;-----------------------------------------------;
+  resmap=True
+;-----------------------------------------------;
+; Get global attributes to set map projection
+;-----------------------------------------------;
+  RPK  = fichier1->RPK ; ZS
+  LON0 = fichier1->LON0
+  BETA = fichier1->BETA
+  if (RPK.gt.0)
+    ;   Lambert  projection from north pole
+    ; ---------------------------
+    resmap@mpProjection          = "LambertConformal"     ; projection
+    resmap@mpLambertParallel1F   =  42  ; to be adjusted !!
+    resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ;  Parallel1 = PArallel2
+    resmap@mpLambertMeridianF    =  LON0      ; get value from globla attribute
+  end if
+
+  if (RPK.lt.0)
+    ;   Lambert projection from south pole
+    ; ---------------------------
+    resmap@mpProjection          = "LambertConformal"     ; projection
+    resmap@mpLambertParallel1F   =  -42  ; to be adjusted !!
+    resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file
+    resmap@mpLambertMeridianF    =  LON0      ; get value from globla attribute
+  end if
+
+  if (RPK.eq.1)
+    ;   Stereographic projection
+    ; ---------------------------
+    resmap@mpProjection = "Stereographic"
+    resmap@mpCenterLonF = LON0
+    resmap@mpCenterRotF = BETA
+    resmap@mpCenterLatF = 90.
+  end if
+
+  if (RPK.eq.-1)
+    ;   Stereographic projection
+    ; ---------------------------
+    resmap@mpProjection = "Stereographic"
+    resmap@mpCenterLonF = LON0
+    resmap@mpCenterRotF = BETA
+    resmap@mpCenterLatF = -90.
+  end if
+
+  if (RPK.eq.0) then
+    ;   Mercator projection
+    ; ---------------------------
+    resmap@mpProjection = "Mercator"
+  end if
+
+ print("Map projection="+resmap@mpProjection)
+
+;====================================
+; Defining the corners for projection
+;====================================
+  resmap@mpLimitMode            = "Corners"
+  resmap@mpLeftCornerLatF       = lat2d(0,0)
+  resmap@mpLeftCornerLonF       = lon2d(0,0)
+  resmap@mpRightCornerLatF     = lat2d(jmax-1,imax-1)
+  resmap@mpRightCornerLonF     = lon2d(jmax-1,imax-1)
+
+;=================================================;
+; lecture des différents champs
+;=================================================;
+  zs  = fichier1->ZS(1:jmax,1:imax) ; ZS
+  zs@lat2d = lat2d
+  zs@lon2d = lon2d
+
+  wt= fichier1->WT(1:kmax+1,1:jmax,1:imax) ; WT
+  vt= fichier1->VT(1:kmax,1:jmax+1,1:imax) ; VT
+  ut= fichier1->UT(1:kmax,1:jmax,1:imax+1) ; UT
+  tht= fichier1->THT(1:kmax,1:jmax,1:imax) ; THT
+  tht@long_name="Potential Temperature"
+  tht@units="K"
+  tht@lat2d = lat2d
+  tht@lon2d = lon2d
+
+  lsthm= fichier1->LSTHM(1:kmax,1:jmax,1:imax) ;LSTHM
+  lsthm@long_name="Large SCale Potential Temperature"
+  lsthm@units="K"
+  lsthm@lat2d = lat2d
+  lsthm@lon2d = lon2d
+
+  lsvm= fichier1->LSVM(1:kmax,1:jmax+1,1:imax) ; LSVM
+  lsvm@long_name="Large SCale Merdian Wind"
+  lsvm@units="m/s"
+  lsvm@lat2d = lat2d
+  lsvm@lon2d = lon2d
+
+  pabst= fichier1->PABST(1:kmax,1:jmax,1:imax) ; PABST
+  pabst@long_name="Pressure"
+  pabst@units="Pa"
+  pabst@lat2d = lat2d
+  pabst@lon2d = lon2d
+
+
+;=================================================;
+; On mets toutes les variables sur la grille 1 
+;=================================================;
+  ut1 = wrf_user_unstagger(ut,"X")
+  ut1@long_name="Zonal wind"
+  ut1@units="m/s"
+  vt1 = wrf_user_unstagger(vt,"Y")
+  vt1@long_name="Meridian wind"
+  vt1@units="m/s"
+  wt1 = wrf_user_unstagger(wt,"Z")
+  wt1@long_name="Vertical wind"
+  wt1@units="m/s"
+
+  ut1@lat2d = lat2d
+  ut1@lon2d = lon2d
+  vt1@lat2d = lat2d
+  vt1@lon2d = lon2d
+  wt1@lat2d = lat2d
+  wt1@lon2d = lon2d
+;=================================================;
+; On calcule l'altitude des champs modèle
+;=================================================;
+
+  zhat= fichier1->ZHAT(1:kmax+1)
+
+; Unstagger zhat (from grid 4 to 1)
+  nzhat=new(kmax,double)
+  do k=0,kmax-1
+    nzhat(k)=(zhat(k)+zhat(k+1))/2.
+  end do
+
+; Create Z3D == ALT
+  alt=new(dimsizes(tht),double)
+  zcoef=1.-zs/nzhat(kmax-1)
+
+  do i=0,imax-1
+    do j=0,jmax-1
+      alt(:,j,i) = nzhat*zcoef(j,i)+zs(j,i)
+    end do
+  end do
+
+;=================================================;
+; Set some other basic resources
+;=================================================;
+  resmap = True
+  resmap@gsnFrame = False
+  resmap@gsnDraw = False
+  resmap@gsnMaximize = True
+  resmap@gsnPaperOrientation = "portrait" 
+  resmap@gsnSpreadColors= True       
+  resmap@tiYAxisString =" "
+  resmap@cnFillOn= True 
+  resmap@cnLinesOn= False 
+
+; If there is an error on HighRes, it means that you don't have the HighRes data
+; You need to download them or the change HighRes by MediumRes
+; See https://www.ncl.ucar.edu/Document/Graphics/rangs.shtml for info
+; -------------------------------------
+  resmap@mpDataBaseVersion     = "HighRes"     ; choose highres map data version (must be donwloaded)
+
+  gsn_define_colormap(wks,"rainbow") 
+
+;=================================================;
+; TRACE
+;=================================================;
+
+; module du vent + vecteur vent pour k=2
+
+; calcul module du vent
+  mod_utvt=sqrt(ut1(1,:,:)*ut1(1,:,:)+vt1(1,:,:)*vt1(1,:,:))
+
+  mod_utvt@long_name="Module of horizontal wind"
+  mod_utvt@lat2d = lat2d
+  mod_utvt@lon2d = lon2d
+  mod_utvt@units="m/s"
+
+; options de tracé
+  res2=resmap
+  res2@cnLevelSelectionMode = "ManualLevels"
+  res2@cnLevelSpacingF    =0.5
+  res2@cnMinLevelValF    = 11
+  res2@cnMaxLevelValF    = 26
+  res2@gsnScalarContour=True
+  res2@vcMinDistanceF = 0.04 
+  res2@vcRefLengthF=0.1      
+  res2@mpFillOn              =  False          ; turn off map fill
+  res2@mpOutlineDrawOrder    = "PostDraw"      ; draw continental outline last
+
+  plot=gsn_csm_vector_scalar_map(wks,ut1(1,:,:),vt1(1,:,:),mod_utvt(:,:),res2)
+  draw(plot)
+  frame(wks)
+
+;  temperature potentielle a 1500m
+  tht_1500m = wrf_user_intrp3d(tht,alt,"h",1500,0.,False)
+  res3=resmap
+  res3@cnLevelSelectionMode = "ManualLevels"
+  res3@cnLevelSpacingF    =0.1
+  res3@cnMinLevelValF    = 302.4
+  res3@cnMaxLevelValF    = 304
+
+  plot_tht = gsn_csm_contour_map(wks,tht_1500m,res3)
+  draw(plot_tht)
+  frame(wks)
+
+
+;  pression +vent  1500m
+  pabst_1500m = wrf_user_intrp3d(pabst,alt,"h",1500,0.,False)
+  ut_1500m = wrf_user_intrp3d(ut1,alt,"h",1500,0.,False)
+  vt_1500m = wrf_user_intrp3d(vt1,alt,"h",1500,0.,False)
+
+  res4=resmap
+  res4@cnLevelSelectionMode = "ManualLevels"
+  res4@cnLevelSpacingF=10
+  res4@cnMinLevelValF= 83800
+  res4@cnMaxLevelValF= 84030
+  res4@gsnScalarContour=True
+  res4@vcMinDistanceF = 0.04 
+  res4@vcRefLengthF=0.1      
+  res4@mpFillOn=  False          
+  res4@mpOutlineDrawOrder= "PostDraw"     
+  res4@vcMinDistanceF = 0.03 
+  res4@vcRefLengthF=0.03     
+
+  plot_pabst=gsn_csm_vector_scalar_map(wks,ut_1500m,vt_1500m,pabst_1500m,res4)
+  draw(plot_pabst)
+  frame(wks)
+
+;=================================================;
+;  Definition  du plan de la coupe verticale
+;=================================================;
+  opts=False
+  plane = new(2,float)
+  plane =(/33,0/)
+  angle=0 
+  X_plane=wrf_user_intrp2d(lat2d,plane,angle,opts)
+  X_desc="latitude"
+
+;=================================================;
+; Coupe verticale de  THT
+;=================================================;
+; inteprolation du champ sur la coupe verticale
+  tht_plane=wrf_user_intrp3d(tht,alt,"v",plane,angle,opts)
+  lsthm_plane=wrf_user_intrp3d(lsthm,alt,"v",plane,angle,opts)
+  vt_plane=wrf_user_intrp3d(vt1,alt,"v",plane,angle,opts)
+  wt_plane=wrf_user_intrp3d(wt1,alt,"v",plane,angle,opts)
+  lsvm_plane=wrf_user_intrp3d(lsvm,alt,"v",plane,angle,opts)
+
+  opts_cv=resmap
+  opts_cv@tiXAxisString           = X_desc
+  opts_cv@tiYAxisString           = "Height (km)"
+  opts_cv@tiXAxisFontHeightF      = 0.020
+  opts_cv@tiYAxisFontHeightF      = 0.020
+  opts_cv@cnMissingValPerimOn     = True
+  opts_cv@cnMissingValFillColor   = "gray"
+  opts_cv@cnMissingValFillPattern = 25
+
+  opts_cv@cnLevelSelectionMode = "ManualLevels"
+  opts_cv@cnLevelSpacingF    =2.5
+  opts_cv@cnMinLevelValF    = 300
+  opts_cv@cnMaxLevelValF    = 355
+  opts_cv@tiMainString="THT"
+  cv_tht = gsn_csm_contour(wks,tht_plane(:,:),opts_cv)
+  draw(cv_tht)
+  frame(wks)
+
+  opts_cv@cnLevelSelectionMode = "ManualLevels"
+  opts_cv@cnLevelSpacingF    =0.1
+  opts_cv@cnMinLevelValF    = -1
+  opts_cv@cnMaxLevelValF    = 1
+  opts_cv@tiMainString="THT-LSTHM"
+  cv_tht_lsthm = gsn_csm_contour(wks,tht_plane(:,:)-lsthm_plane(:,:),opts_cv)
+  draw(cv_tht_lsthm)
+  frame(wks)
+
+  opts_cv@cnLevelSelectionMode = "ManualLevels"
+  opts_cv@cnLevelSpacingF    =0.5
+  opts_cv@cnMinLevelValF    = -5
+  opts_cv@cnMaxLevelValF    = 6
+  opts_cv@tiMainString="VT-LSVM"
+  cv_vt = gsn_csm_contour(wks,vt_plane(:,:)-lsvm_plane(:,:),opts_cv)
+  draw(cv_vt)
+  frame(wks)
+
+  opts_cv@cnLevelSelectionMode = "ManualLevels"
+  opts_cv@cnLevelSpacingF    =0.5
+  opts_cv@cnMinLevelValF    = -6
+  opts_cv@cnMaxLevelValF    = 9
+  opts_cv@tiMainString="WT"
+  cv_wt = gsn_csm_contour(wks,wt_plane(:,:),opts_cv)
+  draw(cv_wt)
+  frame(wks)
+
+end
diff --git a/MY_RUN/KTEST/004_Reunion/006_ncl/run_ncl b/MY_RUN/KTEST/004_Reunion/006_ncl/run_ncl
new file mode 100755 (executable)
index 0000000..74e2d43
--- /dev/null
@@ -0,0 +1,8 @@
+#MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+#MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+#MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+#MNH_LIC for details. version 1.
+set -x
+set -e
+ln -sf ../003_mesonh/REUNI.1.00A20.004.nc4 .
+ncl plot_Reunion.ncl
index 629b0ea..a4e5c76 100644 (file)
@@ -4,7 +4,8 @@ all:
        cd 002_prep_ideal_case && run_prep_ideal_case_xyz 
        cd 003_mesonh          && run_mesonh_xyz 
        cd 004_convdia         && run_conv2dia 
-       cd 005_diaprog         && run_diaprog 
+       cd 005_diaprog         && run_diaprog $
+       cd 006_ncl             && run_ncl
 
 all_ncl_lfi2cdf:
        cd 001_prep_pgd        && get_pgd_files 
@@ -21,6 +22,15 @@ all_ncl_extractdia:
        cd 003_mesonh          && run_mesonh_xyz 
        cd 004_convdia         && run_conv2dia 
        cd 005_ncl_extractdia  && run_ncl 
+
+all_ncl_nc4:
+       cd 001_prep_pgd        && get_pgd_files 
+       cd 001_prep_pgd        && run_prep_pgd_xyz 
+       cd 002_prep_ideal_case && run_prep_ideal_case_xyz 
+       cd 003_mesonh          && run_mesonh_xyz 
+       cd 005_ncl_nc4  && run_ncl 
+       cd 006_ncl             && run_ncl
+
 clean:
        cd 001_prep_pgd        && clean_get_pgd_files 
        cd 001_prep_pgd        && clean_prep_pgd_xyz 
@@ -30,3 +40,5 @@ clean:
        cd 005_diaprog         && clean_diaprog 
        cd 005_ncl_lfi2cdf     && clean_ncl 
        cd 005_ncl_extractdia  && clean_ncl 
+       cd 005_ncl_nc4         && clean_ncl 
+       cd 006_ncl             && clean_ncl