Update mlab script for viz new version. Add stereogaphic proj script.

cmake/SCHISM.local.build

@@ -26,7 +26,7 @@ set( USE_MARSH OFF CACHE BOOLEAN "Use marsh module")
 set( USE_SED2D OFF CACHE BOOLEAN "Use 2D sediment module")
 set( USE_WWM OFF CACHE BOOLEAN "Use wind-wave module")
 set( USE_ICE OFF CACHE BOOLEAN "Use 1-class ICE module")
-set( USE_MICE ON CACHE BOOLEAN "Use multi-class ICE module")
+set( USE_MICE OFF CACHE BOOLEAN "Use multi-class ICE module")
 
 #Tracer models
 set( USE_GEN OFF CACHE BOOLEAN "Use generic tracer module")

src/Utility/Grid_Scripts/arctic_stereo.f90

@@ -0,0 +1,83 @@
+!     Convert between lat/long and Arctic polar stereographic
+!     coordinates for mixed tri/quads.
+!     Input: <fname> (grid format); lat/long must be in degrees
+!     Output: out_<fname>
+!     Depths are not changed.
+!     ifort -O2 -mcmodel=medium -CB -Bstatic -o arctic_stereo arctic_stereo.f90
+      implicit real*8(a-h,o-z)
+      character*70 fname
+      dimension nm(4)
+
+      pi=3.1415926
+      print*, "Input file name:"
+      read(*,'(a70)') fname
+      print*, '1: from lat/long to stereo; -1: from stereo to lat/long:'
+!'
+      read*, ifl
+      if(iabs(ifl).ne.1) then
+        print*, "Incorrect input of ifl",ifl
+        stop
+      endif
+
+      open(12,file=fname,status='old')
+      open(13,file="out_"//fname,status='replace')
+      read(12,*)
+      read(12,*)ne,np
+      write(13,*)
+      write(13,*)ne,np
+      do i=1,np
+        read(12,*)j,xtmp,ytmp,dp
+        call stereographic(ifl,xtmp,ytmp,xout,yout)
+        write(13,'(i10,2(1x,f24.11),1x,e17.8)')j,xout,yout,dp
+      enddo !i
+      do i=1,ne
+        read(12,*)j,k,nm(1:k)
+        write(13,*)j,k,nm(1:k)
+      enddo !i
+      close(12)
+      close(13)
+
+      stop
+      end
+
+!     Polar stereographic projection
+!     LOn/lat inputs/outputs in degrees
+      subroutine stereographic(ifl,xin,yin,xout,yout)
+      implicit real*8(a-h,o-z)
+      parameter(re=6378206.4)
+      parameter(pi=3.1415926)
+      integer, intent(in) :: ifl
+      real*8, intent(in) :: xin,yin
+      real*8, intent(out) :: xout,yout
+
+      if(ifl==1) then !lat/lon to stereo
+        xin1=xin*pi/180 !lon
+        yin1=yin*pi/180
+
+        !3D frame
+        xnd=re*cos(yin1)*cos(xin1)  
+        ynd=re*cos(yin1)*sin(xin1)  
+        znd=re*sin(yin1)
+
+        if(re+znd==0.d0) then
+          print*, 'STEREO: remove south pole,',xin,yin,znd
+          stop
+        endif
+        xout=2*re*xnd/(re+znd) !meters
+        yout=2*re*ynd/(re+znd)
+      else !ifl=-1; stereo to lat/lon
+        rlam=atan2(yin,xin) !radian
+        tmp=2*re*cos(rlam)+xin
+        if(tmp==0.d0) then
+          print*, 'STEREO: south pole encountered,',tmp,xin,yin
+          stop
+        endif
+        half_phi=atan((2*re*cos(rlam)-xin)/tmp)
+
+        xout=rlam/pi*180
+        yout=half_phi*2/pi*180
+      endif
+
+      return
+      end
+

src/Utility/Pre-Processing/list_bnd_pt.f90

@@ -1,4 +1,4 @@
-!     Output open&land bnd points from hgrid.gr3 (for mlab); has
+!     Output open&land&island bnd points from hgrid.gr3 (for mlab); has
 !     duplicate pts btw bnd's
 !     Input: 
 !           (1) hgrid.gr3

src/Utility/Vis_Matlab/get_global_info.m

@@ -1,11 +1,11 @@
 function [ne,np,nvrt,nm,xy00,i23d,ivs,vzcor,vid5,vdry_e,h0,dp,kbp00,vid_eta,nproc,np_lcl,ne_lcl,ns_lcl,iplg,ielg,iegl_rank]=get_global_info(icomb,base,varname,stacks)
 %Get basic info before reading (combined or uncombined) nc outputs 
-%Outputs starting from npoc are only defined if icomb=0 (uncombined)
+%Outputs starting from npoc are only defined if icomb=0 (uncombined).
 %vzcor,vis5,vdry_e: IDs for zcor and varname, and wetdry_elem
 %xy00(np,2); nm(4,ne)
 
 if(icomb==0)
-  fid=fopen([base '/outputs/local_to_global_0000'],'r');
+  fid=fopen([base '/outputs/local_to_global_000000'],'r');
   first=fscanf(fid,'%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n',18);
   %disp(first);
   fclose(fid);
@@ -14,7 +14,7 @@
   nm=NaN(4,ne); i34=zeros(ne,1); xy00=zeros(np,2); iegl_rank=NaN(ne,1);
 
   for irank=0:nproc-1
-    fid=fopen([base '/outputs/local_to_global_' num2str(irank,'%04.f')],'r');
+    fid=fopen([base '/outputs/local_to_global_' num2str(irank,'%06.f')],'r');
     first=fscanf(fid,'%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n',18);
     %char=fscanf(fid,'%s\n')
     ne_lcl(irank+1)=fscanf(fid,'local to global mapping: %d\n',1);
@@ -71,7 +71,7 @@
   end
 
   %Read header of nc
-  fname=[base,'/outputs/' 'schout_0000_' num2str(stacks(1)) '.nc'];
+  fname=[base,'/outputs/' 'schout_000000_' num2str(stacks(1)) '.nc'];
   ncid0 = netcdf.open(fname,'NC_NOWRITE');
   vid5=netcdf.inqVarID(ncid0,varname);
   i23d=netcdf.getAtt(ncid0,vid5,'i23d'); %1: 2D; 2: 3D (whole level)