Gaelle 2/12/2015 : MAJ 009_ICARTT
authorGaelle Tanguy <gaelle.tanguy@meteo.fr>
Wed, 2 Dec 2015 15:39:12 +0000 (15:39 +0000)
committerPhilippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Thu, 19 May 2016 14:44:49 +0000 (16:44 +0200)
16 files changed:
MY_RUN/KTEST/009_ICARTT/001_pgd1/clean_prep_pgd_xyz
MY_RUN/KTEST/009_ICARTT/002_arp2lfi/clean_arp2lfi_xyz
MY_RUN/KTEST/009_ICARTT/002_arp2lfi/run_arp2lfi_xyz
MY_RUN/KTEST/009_ICARTT/003_mesonh/clean_mesonh_xyz
MY_RUN/KTEST/009_ICARTT/004_conv2dia/clean_conv2dia [deleted file]
MY_RUN/KTEST/009_ICARTT/004_conv2dia/dir_conv2dia1 [deleted file]
MY_RUN/KTEST/009_ICARTT/004_conv2dia/run_conv2dia [deleted file]
MY_RUN/KTEST/009_ICARTT/005_cdf/clean_extractdia [deleted file]
MY_RUN/KTEST/009_ICARTT/005_cdf/dir_extractdia1 [deleted file]
MY_RUN/KTEST/009_ICARTT/005_cdf/run_extractdia [deleted file]
MY_RUN/KTEST/009_ICARTT/006_ncl/MESONHtools_ex.ncl [deleted file]
MY_RUN/KTEST/009_ICARTT/006_ncl/clean_ncl
MY_RUN/KTEST/009_ICARTT/006_ncl/plot_ICARTT.ncl [new file with mode: 0644]
MY_RUN/KTEST/009_ICARTT/006_ncl/plot_ICARTT_SEG1.ncl [deleted file]
MY_RUN/KTEST/009_ICARTT/006_ncl/run_ncl
MY_RUN/KTEST/009_ICARTT/Makefile

index 94746d3..0f7d546 100755 (executable)
@@ -1,5 +1,5 @@
 set -x
 
-rm -f *.dir *.hdr *.asc
+rm -f *.dir *.hdr *.asc file_for_xtransfer
 rm -f ICARTT1008_PGD_15km.* OUTPUT_LISTING* OUTPUT_TRANSFER pipe* *.tex
 
index 9aa132e..7f96556 100755 (executable)
@@ -1,5 +1,5 @@
 set -x
 
-rm -f ICART* CPLCH* OUTPUT_LISTING* pipe* 
+rm -f ICART* CPLCH* OUTPUT_LISTING* pipe* file_for_xtransfer 
 rm -f  ecmwf.OD.20040810.18-V2 mocage.GLOB22.20040810.18
 
index 4bf5fb4..47cbdd5 100755 (executable)
@@ -8,6 +8,7 @@ set -e
 
 CHIMIE_FILES=${CHIMIE_FILES:-"$HOME/CHIMIE_FILES"} ; export CHIMIE_FILES
 ln -sf ${CHIMIE_FILES}/*.20040810.18* .
+export GRIB_DEFINITION_PATH=$SRC_MESONH/src/LIB/grib_api-1.13.1/definitions
 
 rm -f ICART* CPLCH20040810.18.* OUTPUT_LISTING* pipe* 
 ln -sf ../001_pgd1/ICARTT1008_PGD_15km.??? .
index a5b725f..fbb097f 100755 (executable)
@@ -1,5 +1,5 @@
 set -x
 rm -f ICART* CPLCH* OUTPUT_LISTING* pipe*
-rm -f EXSEG?.nam PRESSURE 
+rm -f EXSEG?.nam PRESSURE file_for_xtransfer 
 rm -f DATAE1 DATAJ1 DATAS1
 
diff --git a/MY_RUN/KTEST/009_ICARTT/004_conv2dia/clean_conv2dia b/MY_RUN/KTEST/009_ICARTT/004_conv2dia/clean_conv2dia
deleted file mode 100755 (executable)
index ca6739b..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-set -x
-rm -f ICART* 
-rm -f dirconv FICJD LISTING_DIA OUT_DIA pipe_name
diff --git a/MY_RUN/KTEST/009_ICARTT/004_conv2dia/dir_conv2dia1 b/MY_RUN/KTEST/009_ICARTT/004_conv2dia/dir_conv2dia1
deleted file mode 100644 (file)
index e9152f0..0000000
+++ /dev/null
@@ -1,7 +0,0 @@
-2                                                                               
-ICART.1.SEG01.001
-ICART.1.SEG01.002
-ICART.1.SEG01.012cv
-n
-n
-0 
diff --git a/MY_RUN/KTEST/009_ICARTT/004_conv2dia/run_conv2dia b/MY_RUN/KTEST/009_ICARTT/004_conv2dia/run_conv2dia
deleted file mode 100755 (executable)
index ba15f05..0000000
+++ /dev/null
@@ -1,11 +0,0 @@
-#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/ICART.1.SEG01.001.??? .
-ln -sf ../003_mesonh/ICART.1.SEG01.002.??? .
-rm -f ICART.1.SEG01.0*cv.???
-${POSTRUN} conv2dia < dir_conv2dia1
-
diff --git a/MY_RUN/KTEST/009_ICARTT/005_cdf/clean_extractdia b/MY_RUN/KTEST/009_ICARTT/005_cdf/clean_extractdia
deleted file mode 100755 (executable)
index 6107367..0000000
+++ /dev/null
@@ -1,3 +0,0 @@
-set -x
-rm -f ICART* n\'est\ pas\ un\ tty
-rm -f dirextract OUT_DIA
diff --git a/MY_RUN/KTEST/009_ICARTT/005_cdf/dir_extractdia1 b/MY_RUN/KTEST/009_ICARTT/005_cdf/dir_extractdia1
deleted file mode 100644 (file)
index f83e2a5..0000000
+++ /dev/null
@@ -1,62 +0,0 @@
-ICART.1.SEG01.012cv
-ZCDL
-  1
-  0
-  0
-  0
-  0
-  0
-  0
-  0
-  0
-  0
-  0
- 41
- 0.
-250.
-500.
-750.
-1000.
-1250.
-1500.
-1750.
-2000.
-2500.
-2500.
-2750.
-3000.
-3250.
-3500.
-3750.
-4000.
-4250.
-4500.
-4750.
-5000.
-5250.
-5500.
-5750.
-6000.
-6250.
-6500.
-6750.
-7000.
-7350.
-7500.
-7750.
-8000.
-8250.
-8500.
-8750.
-9000.
-9250.
-9500.
-9750.
-10000.
-CONF
-RCT
-RRT
-O3T
-COT
-ZSBIS
-END
diff --git a/MY_RUN/KTEST/009_ICARTT/005_cdf/run_extractdia b/MY_RUN/KTEST/009_ICARTT/005_cdf/run_extractdia
deleted file mode 100755 (executable)
index dca3778..0000000
+++ /dev/null
@@ -1,12 +0,0 @@
-#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 ../004_conv2dia/ICART.1.SEG01.012cv.??? .
-rm -f *.cdf *cdl
-
-extractdia < dir_extractdia1
-
-
diff --git a/MY_RUN/KTEST/009_ICARTT/006_ncl/MESONHtools_ex.ncl b/MY_RUN/KTEST/009_ICARTT/006_ncl/MESONHtools_ex.ncl
deleted file mode 100644 (file)
index 3bec3d5..0000000
+++ /dev/null
@@ -1,262 +0,0 @@
-load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
-
-;*************************************************************************
-; J.-P. CHABOUREAU
-; This is a driver that selects the appropriate 
-; mapping function based upon the file attribute: MAP_PROJ
-; MAP_PROJ=1 [Lambert Conformal]; =2 [Stereographic]; =3 [Mercator]
-;
-; opt: currently not used [potentail use: time counter for XLAT/XLONG]
-;
-; Sample usage:
-;             ncdf = addfile("...", r")
-;             res  = True
-;             WRF_map_c (ncdf, res, 0)
-;             res  = ...
-;
-undef("MESONH_map_c")
-procedure MESONH_map_c (f:file, res:logical, plat, plon, icorner)
-local rank, dimll, nlat, mlon, lat2d, lon2d
-begin
-
-  if (plat(0,0).eq.plat(1,1)) then
-    ZRPK = f->RPK
-;  ZLATOR = f->LATORI
-;  ZLONOR = f->LONORI
-  ZLATOR = f->LATOR
-  ZLONOR = f->LONOR
-;  ZLATOR = 46.59
-;  ZLONOR = -3.09
-  ZBETA = f->BETA
-  ZLAT0 = f->LAT0
-  ZLON0 = f->LON0
-  XHAT = f->W_E_direction
-  YHAT = f->S_N_direction
-  IMAX = dimsizes(XHAT)
-  JMAX = dimsizes(YHAT)
-  zdx=XHAT(2)-XHAT(1)
-  zdy=YHAT(2)-YHAT(1)
- print ("LATOR="+ZLATOR+" - LONOR="+ZLONOR)
- print ("ZLAT0="+ZLAT0+" - ZLON0="+ZLON0)
- print ("ZDX="+zdx+" - RPK="+ZRPK+" - BETA="+ZBETA)
- print ("IMAX="+IMAX+" - JMAX="+JMAX)
-
-;  if (isatt(f,"RPK")) then
-;print ("tutu "+f@RPK)
-;  if (isatt(f,"RPK")) then
-;print ("tutu "+f@RPK)
-;      if (f@MAP_PROJ.eq.1) then
-;          res@mpProjection = "LambertConformal"
-;      end if
-      if (abs(ZRPK).eq.1) then
-        res@mpProjection = "Stereographic"
-        res@mpCenterLonF           = ZLON0
-        res@mpCenterRotF = ZBETA
-        if (ZRPK.eq.1) then
-;;;; XRPK=1: polar stereographic projection from south pole
-          res@mpCenterLatF           = 90. 
-        else
-          res@mpCenterLatF           = -90. 
-        end if
-       end if
-      if (ZRPK.eq.0) then
-;;;; XRPK=0: Mercator projection from earth center
-          res@mpProjection = "Mercator"
-      end if
-      if (abs(ZRPK).gt.0.and.abs(ZRPK).lt.1) then
-;;;; 1>XRPK>0: Lambert projection from south pole
-;;;; -1<XRPK<0: Lambert projection from north pole
-;        res@mpProjection = "LambertConformal"
-        res@mpProjection = "Stereographic"
-        res@mpCenterLonF           = ZLON0
-         if (ZRPK.gt.0) then
-        res@mpCenterLatF           = 90. 
-         else
-        res@mpCenterLatF           = -90. 
-         end if
-      end if
-;  else
-;      print ("MESONH_mapProj: no MAP_PROJ attribute")
-;  end if
-
- print("Map projection="+res@mpProjection)
-
-;=================================================;
-; src/mesonh_MOD/mode_gridproj.f90
-;=================================================;
-  XRADIUS=6371229.0d ; Earth radius (meters)
-  XPI=2.0d*asin(1.)    ; Pi
-  ZRDSDG= XPI/180.0d          ; Radian to Degree conversion factor
-  ZXBM0 = 0.0d
-  ZYBM0 = 0.0d
-
-  lat2d = new((/JMAX,IMAX/),"double")
-  lon2d = new((/JMAX,IMAX/),"double")
-
-;=================================================;
-;  if (abs(ZRPK).eq.1) then
-      if (ZRPK.gt.0.and.ZRPK.le.1) then
-; STEREOGRAPHIC PROJECTION
-    ZCLAT0  = cos(ZRDSDG*ZLAT0)
-    ZSLAT0  = sin(ZRDSDG*ZLAT0)
-    ZCLATOR = cos(ZRDSDG*ZLATOR)
-    ZSLATOR = sin(ZRDSDG*ZLATOR)
-    ZRO0    = (XRADIUS/ZRPK)*(abs(ZCLAT0))^(1.-ZRPK) * \
-              ((1.+ZSLAT0)*abs(ZCLATOR)/(1.+ZSLATOR))^ZRPK
-    ZGA0    = (ZRPK*(ZLONOR-ZLON0)-ZBETA)*ZRDSDG
-    ZXP     = ZXBM0-ZRO0*sin(ZGA0)
-    ZYP     = ZYBM0+ZRO0*cos(ZGA0)
-    do ji=0,IMAX-1
-      do jj=0,JMAX-1
-        ZATA = atan2( -(ZXP-XHAT(ji)) , (ZYP-YHAT(jj)) )/ZRDSDG
-        zlon  = (ZBETA+ZATA)/ZRPK+ZLON0
-        lon2d(jj,ji)=zlon
-        ZRO2 = (XHAT(ji)-ZXP)^2+(YHAT(jj)-ZYP)^2
-        ZJD1 = XRADIUS*(abs(ZCLAT0))^(1.-ZRPK)
-        ZT1  = (ZJD1)^(2./ZRPK)* (1+ZSLAT0)^2
-        ZJD3 = (ZRPK^2*ZRO2)
-        ZT2  = ZJD3
-        ZT2 = ZT2^(1./ZRPK)
-        ZJD1 = (ZT1-ZT2)/(ZT1+ZT2)
-        ZJD1 = acos(ZJD1)
-        ZJD3 = ZJD1
-        zlat = (XPI/2.-ZJD3)/ZRDSDG
-        lat2d(jj,ji)=zlat
-      end do
-    end do
-    
-;   print ("minlat="+min(lat2d)+" maxlat="+max(lat2d))
-;   print ("minlon="+min(lon2d)+" maxlon="+max(lon2d))
-  end if
-;=================================================;
-  if (ZRPK.eq.0) then
-    XBETA=0.
-    XLAT0=0.   ; map reference latitude  (degrees)
-    ZXBM0 = 0.
-    ZYBM0 = 0.
-    ZCGAM    = cos(-ZRDSDG*XBETA)
-    ZSGAM    = sin(-ZRDSDG*XBETA)
-    ZRACLAT0 = XRADIUS*cos(ZRDSDG*ZLAT0)
-    do ji=0,IMAX-1
-      jj=0
-      ZXMI0 = XHAT(ji)-ZXBM0
-      ZYMI0 = YHAT(jj)-ZYBM0
-      zlon  = (ZXMI0*ZCGAM+ZYMI0*ZSGAM)/(ZRACLAT0*ZRDSDG)+ZLONOR
-      do jj=0,JMAX-1
-        lon2d(jj,ji)=zlon
-      end do
-    end do
-    do jj=0,JMAX-1
-      ji=0
-      ZXMI0 = XHAT(ji)-ZXBM0
-      ZYMI0 = YHAT(jj)-ZYBM0
-      ZT1  = log(tan(XPI/4.+ZLATOR*ZRDSDG/2.))
-      ZT2  = (-ZXMI0*ZSGAM+ZYMI0*ZCGAM)/ZRACLAT0
-      zlat  = (-XPI/2.+2.*atan(exp(ZT1+ZT2)))/ZRDSDG
-      do ji=0,IMAX-1
-        lat2d(jj,ji)=zlat
-      end do
-    end do
-    
-  end if
-
-  do ji=0,IMAX-1
-    do jj=0,JMAX-1
-      plat(jj,ji)=lat2d(jj,ji)
-      plon(jj,ji)=lon2d(jj,ji)
-    end do
-  end do
-  end if
-;=================================================;
-  if (icorner(0,0).eq.icorner(1,1)) then
-    icorner(0,0)=0
-    icorner(1,0)=JMAX-1
-    icorner(0,1)=0
-    icorner(1,1)=IMAX-1
-  end if
-; print ("icorner"+icorner)
-;=================================================;
-;  lat2d@units = "degrees_north"    ; not needed
-;  lon2d@units = "degrees_east"
-
-  res@gsnAddCyclic           = False              ; regional data
-
-  res@mpLimitMode            = "Corners"
-  res@mpLeftCornerLatF       = plat(icorner(0,0),icorner(0,1))
-  res@mpLeftCornerLonF       = plon(icorner(0,0),icorner(0,1))
-  res@mpRightCornerLatF     = plat(icorner(1,0),icorner(1,1))
-  res@mpRightCornerLonF     = plon(icorner(1,0),icorner(1,1))
-
-; print ("Corner (0,0); Lat="+res@mpLeftCornerLatF+ \
-;                    ", Lon="+res@mpLeftCornerLonF)
-; print ("Oppos corner; Lat="+res@mpRightCornerLatF+ \
-;                     ", Lon= "+res@mpRightCornerLonF)
-
-;;; print ("0,0 "+XHAT(0)+" lon "+YHAT(0))
-;;; print ("0,0 "+XHAT(IMAX-1)+"  "+YHAT(JMAX-1))
-
-;************************************************
-; Turn on lat / lon labeling
-;************************************************
-  res@pmTickMarkDisplayMode = "Always"         ; turn on tickmarks
-
-  res@mpOutlineBoundarySets = "AllBoundaries" ; state boundaries
-
-  res@mpPerimDrawOrder      = "PostDraw"       ; force map perim
-;************************************************
-; set True for native projection (faster)
-;************************************************
-  res@tfDoNDCOverlay       = True      
-  
-end
-
-;*************************************************************************
-; S. BIELLI
-; This is a routine that interpolate  fields on pressure level for plotting
-; based on pinter.f90 
-; The field to be interpolated must be given at the mass point (grid 1)
-; usage : var_inter=MESONHfunction(var_to_interpol, 850., AbsPressure)
-; Abs pressure must be in Pa
-;
-undef("MESONH_pinter")
-function MESONH_pinter( pfield:numeric, loc_param:numeric, ppabs:numeric )
-
-begin
-
- dimL= dimsizes(loc_param)
-
-; First test for grid = 0
-
-  dimp=dimsizes(ppabs)
-
-  pout=pfield(0:dimL-1,:,:)
-  pfield@_FillValue=999
-  pout@_FillValue=999
-  pout=pout@_FillValue
-
-  do jkp = 0, dimL-1
-       zref=log10(loc_param(jkp)*100.)
-       do jloop = 0, dimp(1)-1
-         do iloop = 0, dimp(2)-1
-            kloop=0
-            flag=True
-           do while (flag .and. (kloop.lt.(dimp(2)-2)))
-             if (.not.ismissing(ppabs(kloop,jloop,iloop))) then
-               zxm=log10(ppabs(kloop,jloop,iloop))
-               zxp=log10(ppabs(kloop+1,jloop,iloop))
-               if ((zxp-zref)*(zref-zxm) .ge. 0) then
-                  pout(jkp,jloop,iloop)= (pfield(kloop,jloop,iloop)*(zxp-zref)+ \
-                       pfield(kloop+1,jloop,iloop)*(zref-zxm))/ (zxp-zxm)
-                  flag=False
-                end if
-              end if
-              kloop=kloop+1
-            end do
-          end do
-        end do
-  end do
-
-  return(pout)
-       
-end
-
index 5603d13..56ddc3f 100755 (executable)
@@ -1,2 +1,2 @@
 set -x
-rm -f zsection_1000.ps
+rm -f zsection_1250.ps *.nc4
diff --git a/MY_RUN/KTEST/009_ICARTT/006_ncl/plot_ICARTT.ncl b/MY_RUN/KTEST/009_ICARTT/006_ncl/plot_ICARTT.ncl
new file mode 100644 (file)
index 0000000..cb4d108
--- /dev/null
@@ -0,0 +1,307 @@
+;================================================;\r
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   \r
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   \r
+; ================================================;\r
+begin\r
+;=================================================;\r
+; open file and read in data\r
+;=================================================;\r
+  a = addfile("ICART.1.SEG01.001.nc4", "r")\r
+  a2 = addfile("ICART.1.SEG01.002.nc4", "r")\r
+\r
+;=================================================;\r
+; Get informations on variable sizes\r
+; dims are dims-2 to remove non-physical values\r
+;=================================================;\r
+  mdims = getfilevardimsizes(a,"THT") ; get dimension sizes\r
+  nd = dimsizes(mdims)\r
+  imax=mdims(nd-1)-2\r
+  jmax=mdims(nd-2)-2\r
+  kmax=mdims(nd-3)-2\r
+\r
+;-------------------------------------------------;\r
+; Read data.\r
+;-------------------------------------------------;\r
+  lat2d = a->LAT(1:jmax,1:imax)\r
+  lat2d@units="degrees_north"\r
+  lon2d = a->LON(1:jmax,1:imax)\r
+  lon2d@units="degrees_east"\r
+\r
+zs  = a->ZS(1:jmax,1:imax) ; ZS\r
+zs@long_name="Orography"\r
+zs@units="m"\r
+zs@lat2d = lat2d\r
+zs@lon2d = lon2d\r
+\r
+printMinMax(zs,0)\r
+\r
+  rc_t1 = a->RCT(1:kmax,1:jmax,1:imax)*1.e3\r
+  rc_t1@long_name="Cloud mixing ratio"\r
+  rc_t1@units="g/kg"\r
+  rc_t1@lat2d=lat2d\r
+  rc_t1@lon2d=lon2d\r
+printMinMax(rc_t1,0)\r
+\r
+;\r
+  o3_t1 = a->O3T(1:kmax,1:jmax,1:imax)*1.e9\r
+  o3_t1@long_name="Ozone"\r
+  o3_t1@units="ppbv"\r
+  o3_t1@lat2d=lat2d\r
+  o3_t1@lon2d=lon2d\r
+\r
+;\r
+  co_t1 = a->COT(1:kmax,1:jmax,1:imax)*1.e9\r
+  co_t1@long_name="carbon monoxide"\r
+  co_t1@units="ppbv"\r
+  co_t1@lat2d=lat2d\r
+  co_t1@lon2d=lon2d\r
+;\r
+;\r
+  rc_t2 = a2->RCT(1:kmax,1:jmax,1:imax)*1.e3\r
+  rc_t2@long_name="Cloud mixing ratio"\r
+  rc_t2@units="g/kg"\r
+  rc_t2@lat2d=lat2d\r
+  rc_t2@lon2d=lon2d\r
+\r
+;\r
+  o3_t2 = a2->O3T(1:kmax,1:jmax,1:imax)*1.e9\r
+  o3_t2@long_name="Ozone"\r
+  o3_t2@units="ppbv"\r
+  o3_t2@lat2d=lat2d\r
+  o3_t2@lon2d=lon2d\r
+\r
+;\r
+  co_t2 = a2->COT(1:kmax,1:jmax,1:imax)*1.e9\r
+  co_t2@long_name="carbon monoxide"\r
+  co_t2@units="ppbv"\r
+  co_t2@lat2d=lat2d\r
+  co_t2@lon2d=lon2d\r
+\r
+;-----------------------------------------------;\r
+;=================================================;\r
+; On calcule l'altitude des champs modèle\r
+;=================================================;\r
+\r
+zhat= a->ZHAT(1:kmax+1)\r
+\r
+; Unstagger zhat (from grid 4 to 1)\r
+    nzhat=new(kmax,double)\r
+    do k=0,kmax-1\r
+     nzhat(k)=(zhat(k)+zhat(k+1))/2.\r
+    end do\r
+\r
+; Create Z3D == ALT\r
+    alt=new(dimsizes(rc_t1),double)\r
+    zcoef=1.-zs/nzhat(kmax-1)\r
+\r
+    do i=0,imax-1\r
+    do j=0,jmax-1\r
+       alt(:,j,i) = nzhat*zcoef(j,i)+zs(j,i)\r
+    end do\r
+    end do\r
+\r
+alt@lat2d = lat2d\r
+alt@lon2d = lon2d\r
+\r
+\r
+\r
+\r
+;-----------------------------------------------;\r
+; Set map projection ressources using projection parameters\r
+;-----------------------------------------------;\r
+; Read projection parameters\r
+; --------------------\r
+    RPK  = a->RPK\r
+    BETA = a->BETA\r
+    LON0 = a->LON0\r
+\r
+  resmap=True\r
+  if (RPK.gt.0)\r
+; ---------------------------\r
+  ;   Lambert  projection from north pole\r
+; ---------------------------\r
+   resmap@mpProjection          = "LambertConformal"     ; projection\r
+   pole                         = 1    ; projection for north hemisphere\r
+   resmap@mpLambertParallel1F   = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere\r
+   resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file\r
+   resmap@mpLambertMeridianF    = LON0      ; ncl adds from grib file\r
+  end if\r
+\r
+  if (RPK.lt.0)\r
+; ---------------------------\r
+  ;   Lambert projection from south pole\r
+; ---------------------------\r
+   resmap@mpProjection          = "LambertConformal"     ; projection\r
+   pole                         = -1                     ; projection for south hemisphere\r
+   resmap@mpLambertParallel1F   = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere\r
+   resmap@mpLambertParallel2F   = resmap@mpLambertParallel1F  ; ncl adds from grib file\r
+   resmap@mpLambertMeridianF    =  LON0      ; ncl adds from grib file\r
+  end if\r
+\r
+  if (RPK.eq.1)\r
+; ---------------------------\r
+  ;   Stereographic projection north\r
+; ---------------------------\r
+    resmap@mpProjection = "Stereographic"\r
+    resmap@mpCenterLonF           = LON0\r
+    resmap@mpCenterRotF           = BETA\r
+    resmap@mpCenterLatF           = 90\r
+  end if\r
+\r
+  if (RPK.eq.-1)\r
+; ---------------------------\r
+  ;   Stereographic projection south\r
+; ---------------------------\r
+    resmap@mpProjection = "Stereographic"\r
+    resmap@mpCenterLonF           = LON0\r
+    resmap@mpCenterRotF           = BETA\r
+    resmap@mpCenterLatF           = -90\r
+  end if\r
+\r
+  if (RPK.eq.0) then\r
+; ---------------------------\r
+  ;   Mercator projection\r
+; ---------------------------\r
+    resmap@mpProjection = "Mercator"\r
+  end if\r
+\r
+ print("Map projection="+resmap@mpProjection)\r
+\r
+; Defining the corners for projection\r
+; --------------------------------\r
+  resmap@mpLimitMode            = "Corners"\r
+  resmap@mpLeftCornerLatF       = lat2d(0,0)\r
+  resmap@mpLeftCornerLonF       = lon2d(0,0)\r
+  resmap@mpRightCornerLatF     = lat2d(jmax-1,imax-1)\r
+  resmap@mpRightCornerLonF     = lon2d(jmax-1,imax-1)\r
+  \r
+;=================================================;\r
+; PLOT\r
+;=================================================;\r
+; interpolation des champs a 1250 m\r
+rc_t1_plane = wrf_user_intrp3d(rc_t1,alt,"h",1250,0.,False)\r
+printMinMax(rc_t1_plane,0)\r
+printMinMax(alt,0)\r
+\r
+rc_t2_plane = wrf_user_intrp3d(rc_t2,alt,"h",1250,0.,False)\r
+co_t1_plane = wrf_user_intrp3d(co_t1,alt,"h",1250,0.,False)\r
+co_t2_plane = wrf_user_intrp3d(co_t2,alt,"h",1250,0.,False)\r
+o3_t1_plane = wrf_user_intrp3d(o3_t1,alt,"h",1250,0.,False)\r
+o3_t2_plane = wrf_user_intrp3d(o3_t2,alt,"h",1250,0.,False)\r
+\r
+\r
+  figname ="zsection_1250"\r
+  wks  = gsn_open_wks("ps",figname)   ; open a ncgm file\r
+  gsn_define_colormap(wks,"WhBlGrYeRe") ; Choose colormap\r
+\r
+  res                 = resmap          \r
+  res@gsnDraw                  = False         ; don't draw yet\r
+  res@gsnFrame                 = False         ; don't advance frame yet\r
+\r
+; X-axis title (tiY)                              \r
+  res@tiXAxisFontHeightF = 0.018  ; font height\r
+  res@tiXAxisFont        = 21     ; font index\r
+  res@tiXAxisString      = "longitude"  ; string to use as the X-Axis title\r
+\r
+; Y-axis title (tiY)\r
+  res@tiYAxisFontHeightF = 0.018  ; font height\r
+  res@tiYAxisFont        = 21     ; font index\r
+  res@tiYAxisString      = "latitude"  ; string to use as the Y-Axis title\r
+\r
+; BW\r
+  res@cnLinesOn            = False           \r
+  res@cnFillOn             = True        \r
+  res@gsnSpreadColors      = True \r
+;\r
+; label bar (lb)\r
+;  res@lbAutoManage       = False  \r
+;  res@lbBottomMarginF    = 0.4        ; offset\r
+;  res@lbOrientation      = "Vertical"\r
+\r
+; Map ressources \r
+;  res@mpDataBaseVersion       = "HighRes"     ; choose highres map data version (must be donwloaded)\r
+;  res@mpDataBaseVersion       = "MediumRes"   ; choose highres map data version (must be donwloaded)\r
+  res@mpGridAndLimbOn          = True             ; turn on lat/lon lines\r
+  res@mpGridLatSpacingF        = 10.              ; spacing for lat lines\r
+  res@mpGridLonSpacingF        = 10.              ; spacing for lon lines\r
+\r
+  res@mpGeophysicalLineColor = "Black"  ; default value in lowres\r
+  res@mpNationalLineColor    = "Black"  ; idem\r
+  res@mpUSStateLineColor     = "Black"  ; idem\r
+  res@mpGridLineColor        = "Black"\r
+  res@mpLimbLineColor        = "Black"\r
+  res@mpPerimLineColor       = "Black"\r
+\r
+  \r
+  res@gsnCenterString="heure=19"\r
+\r
+; plot cloud mixing ratio\r
+  res@cnLevelSelectionMode = "ExplicitLevels"\r
+  res@cnLevels =     (/0.01,0.015,0.02,0.025,0.03,0.035,0.04,0.045,0.05,0.055,0.06/)\r
+  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,96,101/) ; color of a contour\r
+  plot_rc = gsn_csm_contour_map(wks,rc_t1_plane(:,:),res)\r
+  draw(plot_rc)\r
+  frame(wks)\r
+  delete(res@cnLevels)\r
+  delete(res@cnFillColors)\r
+\r
+; plot ozone\r
+  res@cnLevelSelectionMode = "ExplicitLevels" ; method for selecting the contour\r
+  res@cnLevels     = (/15.,20., 25., 35., 40., 45., 50., 55., 60., 65./)\r
+  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,96/) ; color of a contour\r
+;  res@cnLevelSelectionMode = "AutomaticLevels" \r
+  plot_o3 = gsn_csm_contour_map(wks,o3_t1_plane(:,:),res)\r
+  draw(plot_o3)\r
+  frame(wks)\r
+  delete(res@cnLevels)\r
+  delete(res@cnFillColors)\r
+\r
+; plot co\r
+  res@cnLevelSelectionMode = "ExplicitLevels" ; method for selecting the contour\r
+  res@cnLevels     = (/110.,112.5,115.,117.5,120.,122.5,125.,127.5,130.,132.5,135./)\r
+  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,96,101/) ; color of a contour\r
+;  res@cnLevelSelectionMode = "AutomaticLevels" \r
+  plot_co = gsn_csm_contour_map(wks,co_t1_plane(:,:),res)\r
+  draw(plot_co)\r
+  frame(wks)\r
+  delete(res@cnLevels)\r
+  delete(res@cnFillColors)\r
+\r
+  res@gsnCenterString="heure=20"\r
+\r
+; plot cloud mixing ratio\r
+  res@cnLevelSelectionMode = "ExplicitLevels"\r
+  res@cnLevels =     (/0.01,0.015,0.02,0.025,0.03,0.035,0.04,0.045,0.05,0.055,0.06/)\r
+  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,96,101/) ; color of a contour\r
+  plot_rc1 = gsn_csm_contour_map(wks,rc_t2_plane(:,:),res)\r
+  draw(plot_rc1)\r
+  frame(wks)\r
+  delete(res@cnLevels)\r
+  delete(res@cnFillColors)\r
+\r
+; plot ozone\r
+  res@cnLevelSelectionMode = "ExplicitLevels" ; method for selecting the contour\r
+  res@cnLevels     = (/15.,20., 25., 35., 40., 45., 50., 55., 60., 65./)\r
+  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,96/) ; color of a contour\r
+;  res@cnLevelSelectionMode = "AutomaticLevels" \r
+  plot_o31 = gsn_csm_contour_map(wks,o3_t2_plane(:,:),res)\r
+  draw(plot_o31)\r
+  frame(wks)\r
+  delete(res@cnLevels)\r
+  delete(res@cnFillColors)\r
+\r
+; plot co\r
+  res@cnLevelSelectionMode = "ExplicitLevels" ; method for selecting the contour\r
+  res@cnLevels     = (/110.,112.5,115.,117.5,120.,122.5,125.,127.5,130.,132.5,135./)\r
+  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,96,101/) ; color of a contour\r
+;  res@cnLevelSelectionMode = "AutomaticLevels" \r
+  plot_co1 = gsn_csm_contour_map(wks,co_t2_plane(:,:),res)\r
+  draw(plot_co1)\r
+  frame(wks)\r
+\r
+;;;;;;;;;;;;;;;;;;;;;;;;\r
+\r
+end\r
+\r
+\r
+\r
diff --git a/MY_RUN/KTEST/009_ICARTT/006_ncl/plot_ICARTT_SEG1.ncl b/MY_RUN/KTEST/009_ICARTT/006_ncl/plot_ICARTT_SEG1.ncl
deleted file mode 100644 (file)
index 55a4885..0000000
+++ /dev/null
@@ -1,181 +0,0 @@
-;================================================;
-load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
-load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
-load "MESONHtools_ex.ncl"   
-; ================================================;
-begin
-;=================================================;
-; open file and read in data
-;=================================================;
-  dir_input1="$SRC_MESONH/MY_RUN/KTEST/009_ICARTT/005_cdf"
-  a = addfile(dir_input1+"/ICART.1.SEG01.012cvZCL.nc", "r")
-;-------------------------------------------------;
-; Time values and units.
-;-------------------------------------------------;
-  time = a ->time
-  time@units = "seconds since 2002-01-0100:00:00"
-; Convert to UTC time.
-  utc_date = ut_calendar(time, 0)
-; Store return information into more meaningful variables.
-  year   = floattointeger(utc_date(:,0))    ; Convert to integer for
-  month  = floattointeger(utc_date(:,1))    ; use sprinti
-  day    = floattointeger(utc_date(:,2))
-  hour   = floattointeger(utc_date(:,3))
-  minute = floattointeger(utc_date(:,4))
-  second = utc_date(:,5)
-; Array to hold month abbreviations. Don't store anything in index
-; '0' (i.e. let index 1=Jan, 2=Feb, ..., index 12=Dec).
-  month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
-                   "Oct","Nov","Dec"/)
-  date_str = sprinti("%0.2iZ ", hour) + sprinti("%0.2i ", day) + \
-             month_abbr(month) + " "  + sprinti("%0.4i", year)
-  hh = sprinti("%0.2i", hour)
-;-------------------------------------------------;
-; Read data.
-;-------------------------------------------------;
-  rc = a->RCT(:,:,:,:)*1.e3
-  rc@long_name="Cloud mixing ratio"
-  rc@units="g/kg"
-  rc!0 = "time"
-  rc&time=hh
-  rc!1 = "alt"
-  rc&alt = a->vertical_levels(:)
-;
-  o3 = a->O3T(:,:,:,:)*1.e9
-  o3@long_name="Ozone"
-  o3@units="ppbv"
-  o3!0="time"
-  o3&time=hh
-  o3!1 = "alt"
-  o3&alt = a->vertical_levels(:)
-;
-  co = a->COT(:,:,:,:)*1.e9
-  co@long_name="carbon monoxide"
-  co@units="ppbv"
-  co!0="time"
-  co&time=hh
-  co!1 = "alt"
-  co&alt = a->vertical_levels(:)
-;
-;  nlevel = 7 ; 1500 m
-  nlevel = 5  ; 1000 m
-  
-;=================================================;
-; PLOT
-;=================================================;
-
-  figname ="zsection_1000"
-  wks  = gsn_open_wks("ps",figname)   ; open a ncgm file
-  gsn_define_colormap(wks,"WhBlGrYeRe") ; Choose colormap
-
-  res                 = True          
-  res@gsnDraw                  = False         ; don't draw yet
-  res@gsnFrame                 = False         ; don't advance frame yet
-
-; X-axis title (tiY)                              
-  res@tiXAxisFontHeightF = 0.018  ; font height
-  res@tiXAxisFont        = 21     ; font index
-  res@tiXAxisString      = "longitude"  ; string to use as the X-Axis title
-
-; Y-axis title (tiY)
-  res@tiYAxisFontHeightF = 0.018  ; font height
-  res@tiYAxisFont        = 21     ; font index
-  res@tiYAxisString      = "latitude"  ; string to use as the Y-Axis title
-
-; BW
-  res@cnLinesOn            = False           
-  res@cnFillOn             = True        
-  res@gsnSpreadColors      = True 
-;
-; label bar (lb)
-;  res@lbAutoManage       = False  
-;  res@lbBottomMarginF    = 0.4        ; offset
-;  res@lbOrientation      = "Vertical"
-
-; Map ressources 
-;  res@mpDataBaseVersion       = "HighRes"     ; choose highres map data version (must be donwloaded)
-;  res@mpDataBaseVersion       = "MediumRes"   ; choose highres map data version (must be donwloaded)
-  res@mpGridAndLimbOn          = True             ; turn on lat/lon lines
-  res@mpGridLatSpacingF        = 10.              ; spacing for lat lines
-  res@mpGridLonSpacingF        = 10.              ; spacing for lon lines
-
-  res@mpGeophysicalLineColor = "Black"  ; default value in lowres
-  res@mpNationalLineColor    = "Black"  ; idem
-  res@mpUSStateLineColor     = "Black"  ; idem
-  res@mpGridLineColor        = "Black"
-  res@mpLimbLineColor        = "Black"
-  res@mpPerimLineColor       = "Black"
-
-; Use MESONHtools procedure to set map resources
-  dimrc=dimsizes(rc)
-  IMAX=dimrc(3)
-  JMAX=dimrc(2)  
-  latmod1 = new((/JMAX,IMAX/),"double")
-  latmod1(:,:)=0.
-  lonmod1 = new((/JMAX,IMAX/),"double")
-  cormod1 = new((/2,2/),"integer")
-  cormod1(:,:)=0
-  MESONH_map_c (a, res, latmod1, lonmod1, cormod1)  
-  
-  res@gsnCenterString="heure="+hh(0)
-
-; plot cloud mixing ratio
-  res@cnLevelSelectionMode = "AutomaticLevels" 
-;  res@cnLevelSelectionMode = "ExplicitLevels"
-;  res@cnLevels =     (/1.e-5/)
-  plot_rc = gsn_csm_contour_map(wks,rc(0,nlevel,:,:),res)
-  draw(plot_rc)
-  frame(wks)
-
-; plot ozone
-  res@cnLevelSelectionMode = "ExplicitLevels" ; method for selecting the contour
-  res@cnLevels     = (/15.,20., 25., 35., 40., 45., 50., 55., 60., 65./)
-  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,102,2/) ; color of a contour
-;  res@cnLevelSelectionMode = "AutomaticLevels" 
-  plot_o3 = gsn_csm_contour_map(wks,o3(0,nlevel,:,:),res)
-  draw(plot_o3)
-  frame(wks)
-
-; plot co
-;  res@cnLevelSelectionMode = "ExplicitLevels" ; method for selecting the contour
-;  res@cnLevels     = (/15.,20., 25., 35., 40., 45., 50., 55., 60., 65./)
-;  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,102,2/) ; color of a contour
-  res@cnLevelSelectionMode = "AutomaticLevels" 
-  plot_co = gsn_csm_contour_map(wks,co(0,nlevel,:,:),res)
-  draw(plot_co)
-  frame(wks)
-
-  res@gsnCenterString="heure="+hh(1)
-
-; plot cloud mixing ratio
-  res@cnLevelSelectionMode = "AutomaticLevels" 
-;  res@cnLevelSelectionMode = "ExplicitLevels"
-;  res@cnLevels =     (/1.e-5/)
-  plot_rc1 = gsn_csm_contour_map(wks,rc(1,nlevel,:,:),res)
-  draw(plot_rc1)
-  frame(wks)
-
-; plot ozone
-  res@cnLevelSelectionMode = "ExplicitLevels" ; method for selecting the contour
-  res@cnLevels     = (/15.,20., 25., 35., 40., 45., 50., 55., 60., 65./)
-  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,102,2/) ; color of a contour
-;  res@cnLevelSelectionMode = "AutomaticLevels" 
-  plot_o31 = gsn_csm_contour_map(wks,o3(1,nlevel,:,:),res)
-  draw(plot_o31)
-  frame(wks)
-
-; plot co
-;  res@cnLevelSelectionMode = "ExplicitLevels" ; method for selecting the contour
-;  res@cnLevels     = (/15.,20., 25., 35., 40., 45., 50., 55., 60., 65./)
-;  res@cnFillColors = (/2,6,12,40,45,51,62,72,80,89,102,2/) ; color of a contour
-  res@cnLevelSelectionMode = "AutomaticLevels" 
-  plot_co1 = gsn_csm_contour_map(wks,co(1,nlevel,:,:),res)
-  draw(plot_co1)
-  frame(wks)
-
-;;;;;;;;;;;;;;;;;;;;;;;;
-
-end
-
-
-
index a6136d0..09f1e79 100755 (executable)
@@ -5,6 +5,8 @@
 set -x
 set -e
 
+ln -sf ../003_mesonh/ICART.1.SEG01.001.nc4 .
+ln -sf ../003_mesonh/ICART.1.SEG01.002.nc4 .
+
 rm -f *.ps 
-ncl plot_ICARTT_SEG1.ncl
-gs zsection_1000.ps
+ncl plot_ICARTT.ncl
index 16b39a2..53eafb8 100644 (file)
@@ -1,4 +1,4 @@
-all: E001_pgd1 E002_arp2lfi E003_mesonh E004_conv2dia E005_cdf E006_ncl 
+all: E001_pgd1 E002_arp2lfi E003_mesonh  E006_ncl 
 
 E001_pgd1 :
        cd 001_pgd1            && get_chimie_files
@@ -7,10 +7,6 @@ E002_arp2lfi :
        cd 002_arp2lfi         && run_arp2lfi_xyz
 E003_mesonh:
        cd 003_mesonh          && run_mesonh_xyz
-E004_conv2dia:
-       cd 004_conv2dia        && run_conv2dia
-E005_cdf:
-       cd 005_cdf             && run_extractdia
 E006_ncl:
        cd 006_ncl             && run_ncl
 
@@ -18,7 +14,5 @@ clean:
        cd 001_pgd1            && clean_prep_pgd_xyz
        cd 002_arp2lfi         && clean_arp2lfi_xyz
        cd 003_mesonh          && clean_mesonh_xyz
-       cd 004_conv2dia        && clean_conv2dia
-       cd 005_cdf             && clean_extractdia
        cd 006_ncl             && clean_ncl