;; Example script for creating a file o rotation angles to be used in CliMAF with ;; operator 'plot' (i.e. with script gplot.ncl) ;; This one iwas created by Matthieu Chevallier, and has been tested with ORCA1 grid ;Step 1. Load the necessary NCL scripts 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" ; ; ; FONCTION AUXILIAIRE. ; function cl_recouvrement(field,psg) local ntot,jpj2m,ll1,ll2,dimse,dimll1,dimll2,l1tmp,l2tmp,l3tmp,l4tmp begin dimse=dimsizes(field) ntot=dimse(2) jpj2m=ntot/2 ll1=field(0,dimse(1)-1-1,:) ll2=field(0,dimse(1)-1-1,:) dimll1=dimsizes(ll1) dimll2=dimsizes(ll2) l2tmp=ll2(jpj2m-1:dimll2-1-1) l1tmp=ll1(jpj2m-1:dimll1-1-1) l3tmp=ll2(0:jpj2m-1) l4tmp=ll1(0:jpj2m-1) field(0,dimse(1)-1-1,:jpj2m-1)=psg*l2tmp(::-1) field(0,dimse(1)-2-1,:jpj2m-1)=psg*l1tmp(::-1) field(0,dimse(1)-1-1,jpj2m:)=psg*l3tmp(::-1) field(0,dimse(1)-2-1,jpj2m:)=psg*l4tmp(::-1) return(field) ;# end ; ; =========================================================== ; begin ; ======================================================== conf="ORCA1" if conf.eq."ORCA12" then repfile="/dataref/rd/INITIALISATION/ORCA12/T321/ORCA12-T321.mesh_hgr.nc" nx=3059 ny=4322 else if conf.eq."ORCA025" then repfile="/dataref/rd/INITIALISATION/ORCA025_LIM/ORCA025_LIM-T322/ORCA025_LIM-T322_mesh_hgr.nc" nx=1021 ny=1442 end if end if if conf.eq."ORCA1" then repfile="/cnrm/est/COMMON/CDAT5/Data/mesh_mask_NEMO1.2.nc" nx=292 ny=362 end if ; ;fmesh= addfile("/sxaster1/data1/UTILS/CDAT5/Data/mesh_mask_NEMO1.2.nc","r") fmesh= addfile(repfile,"r") glamf = fmesh->glamf ;printVarSummary(glamf) ;ax=glamf&x ;at=glamf&t ;ay=glamf&y ;print(top) gphif = fmesh->gphif glamu = fmesh->glamu glamv = fmesh->glamv gphiu = fmesh->gphiu gphiv = fmesh->gphiv ; print("Champs charges... en avant!") ; ; Les champs sont sur la grille du modele. ; On les tourne. ; Ces operations sont extraites de rot_rep (NEMO). ; pi=4*atan(1) rad=pi/180. ; ;(xnpu,ynpu,nnpu) = np_dir_mod(glamu,gphiu) ; glamu0=glamu(0,:,:) gphiu0=gphiu(0,:,:) xnpu = 0. - 2.*cos(rad*glamu0)*tan(pi/4 - rad*gphiu0/2) ynpu = 0. - 2.*sin(rad*glamu0)*tan(pi/4 - rad*gphiu0/2) nnpu = xnpu*xnpu + ynpu*ynpu ; ;(xnpv,ynpv,nnpv) = np_dir_mod(glamv,gphiv) ; glamv0=glamv(0,:,:) gphiv0=gphiv(0,:,:) xnpv = 0. - 2.*cos(rad*glamv0)*tan(pi/4 - rad*gphiv0/2) ynpv = 0. - 2.*sin(rad*glamv0)*tan(pi/4 - rad*gphiv0/2) nnpv = xnpv*xnpv + ynpv*ynpv ; ;(xffu,yffu,nffu) = dir_func(glamf,gphif,nnpu,0,1) ; ; ;delete(dime) dime=dimsizes(glamf) xffu=new((/dime(1),dime(2)/),float,-9999.) yffu=new((/dime(1),dime(2)/),float,-9999.) nffu=new((/dime(1),dime(2)/),float,-9999.) glam=glamf(0,1:, :) glan=glamf(0,0:dime(1)-1-1,:) gphi=gphif(0,1:, :) gphh=gphif(0,0:dime(1)-1-1,:) xffu(1:,:) = 2.*cos(rad*glam)*tan(pi/4 - rad*gphi/2) \ - 2.*cos(rad*glan)*tan(pi/4 - rad*gphh/2) yffu(1:,:) = 2.*sin(rad*glam)*tan(pi/4 - rad*gphi/2) \ - 2.*sin(rad*glan)*tan(pi/4 - rad*gphh/2) nffu = sqrt(nnpu * (xffu*xffu + yffu*yffu)) nffu=where(nffu.le.1e-14,1e-14,nffu) ; ; (xffv,yffv,nffv) = dir_func(glamf,gphif,nnpv,1,0) ; delete(glam) delete(glan) delete(gphi) delete(gphh) glam=glamf(0,:, 1:) glan=glamf(0,:,0:dime(2)-1-1) gphi=gphif(0,:, 1:) gphh=gphif(0,:,0:dime(2)-1-1) ; xffv=new((/dime(1),dime(2)/),float,-9999.) yffv=new((/dime(1),dime(2)/),float,-9999.) nffv=new((/dime(1),dime(2)/),float,-9999.) xffv(:,1:) = 2.*cos(rad*glam)*tan(pi/4 - rad*gphi/2) \ - 2.*cos(rad*glan)*tan(pi/4 - rad*gphh/2) yffv(:,1:) = 2.*sin(rad*glam)*tan(pi/4 - rad*gphi/2) \ - 2.*sin(rad*glan)*tan(pi/4 - rad*gphh/2) ; nffv = sqrt(nnpv * (xffv*xffv + yffv*yffv)) nffv=where(nffv.le.1e-14,1e-14,nffv) ; ;# cosinus and sinus using scalar and vectorial products gsinu = ( xnpu*yffu - ynpu*xffu ) / nffu gcosu = ( xnpu*xffu + ynpu*yffu ) / nffu gsinv = ( xnpv*xffv + ynpv*yffv ) / nffv gcosv = - ( xnpv*yffv - ynpv*xffv ) / nffv ; ;#GEOGRAPHICAL MESH ;#COSU - SINU ; (gsinu,gcosu,gdgla) = geo_mesh(glamf,gsinu,gcosu,0,1) delete(dime) delete(glam) delete(glan) dime=dimsizes(glamf) dgla = new((/dime(1),dime(2)/),float,-9999.) glam = glamf(0,1:, :) glan = glamf(0,0:dime(1)-1-1,:) dgla(1:,:) = mod(abs(glam-glan),360) gsinu = where(dgla.lt.1e-3,0., gsinu) gcosu = where(dgla.lt.1e-3,1., gcosv) ; ;gsinu = mask(gsinu,abs(gsinu).gt.2,False) ;gcosu = mask(gcosu,abs(gcosu).gt.2,False) ;gsinu = where(gsinu.ge.1,1.,gsinu) ;gcosu = where(gcosu.ge.1,1.,gcosu) ;gsinu = where(gsinu.le.-1,-1.,gsinu) ;gcosu = where(gcosu.le.-1,-1.,gcosu) gsinu@_FillValue = -9999. gcosu@_FillValue = -9999. ; gsinur=reshape(gsinu,(/1,nx,ny/)) gcosur=reshape(gcosu,(/1,nx,ny/)) ; ;gsinur=cl_recouvrement(gsinur,1) ;gcosur=cl_recouvrement(gcosur,1) print("OK COSU et SINU!") ;#COSV - SINV ; dimphif=dimsizes(gphif) gphim = gphif(0,:,1:) gphin = gphif(0,:,0:dimphif(2)-1-1) gladiffmod=new((/dimphif(1),dimphif(2)/),float,-9999.) gladiffmod(:,1:)=abs(gphim - gphin) ;#gladiffmod = abs(gphim - gphin) gsinv = where(gladiffmod.lt.1e-3,0., gsinv) gcosv = where(gladiffmod.lt.1e-3,1., gcosv) ; ;gsinv = mask(gsinv,abs(gsinv).gt.2,False) ;gcosv = mask(gcosv,abs(gcosv).gt.2,False) ;gsinv = where(gsinv.ge.1,1.,gsinv) ;gcosv = where(gcosv.ge.1,1.,gcosv) ;gsinv = where(gsinv.le.-1,-1.,gsinv) ;gcosv = where(gcosv.le.-1,-1.,gcosv) gsinv@_FillValue = -9999. gcosv@_FillValue = -9999. ; gsinvr=reshape(gsinv,(/1,nx,ny/)) gcosvr=reshape(gcosv,(/1,nx,ny/)) ;gsinvr=cl_recouvrement(gsinvr,1) ;gcosvr=cl_recouvrement(gcosvr,1) ; gsinvr@short_name="GSINV" gcosvr@short_name="GCOSV" gsinur@short_name="GSINU" gcosur@short_name="GCOSU" ; gsinvr@long_name="GSINV" gcosvr@long_name="GCOSV" gsinur@long_name="GSINU" gcosur@long_name="GCOSU" ; gsinvr@missing_value=-9999. gcosvr@missing_value=-9999. gsinur@missing_value=-9999. gcosur@missing_value=-9999. ; gsinvr@valid_min=-1. gcosvr@valid_min=-1. gsinur@valid_min=-1. gcosur@valid_min=-1. ; gsinvr@valid_max=1. gcosvr@valid_max=1. gsinur@valid_max=1. gcosur@valid_max=1. ; gsinvr@axis="TYX" gcosvr@axis="TYX" gsinur@axis="TYX" gcosur@axis="TYX" printVarSummary(gcosur) printVarSummary(glamf) print("OK COSV et SINV!") ; ======================================================== print("ON ECRIT.") fon="angle_"+conf+".nc" fo = addfile(fon,"c") setfileoption(fo,"DefineMode",True) dimNames = (/"t","y", "x"/) dimSizes = (/1, nx, ny/) dimUnlim=(/True,False,False/) filedimdef(fo,dimNames,dimSizes,dimUnlim) print("Var GCOSU") filevardef(fo,"GCOSU",typeof(glamf),getvardims(glamf)) print("Att GCOSU") filevarattdef(fo,"GCOSU",gcosur) print("Wri GCOSU") fo->GCOSU = (/gcosur/) filevardef(fo,"GCOSV",typeof(glamf),getvardims(glamf)) filevarattdef(fo,"GCOSV",gcosvr) fo->GCOSV = (/gcosvr/) filevardef(fo,"GSINU",typeof(glamf),getvardims(glamf)) filevarattdef(fo,"GSINU",gsinur) fo->GSINU = (/gsinur/) filevardef(fo,"GSINV",typeof(glamf),getvardims(glamf)) filevarattdef(fo,"GSINV",gsinvr) fo->GSINV = (/gsinvr/) ; end