/[MITgcm]/MITgcm_contrib/mlosch/interp_llc/addshelfice.m
ViewVC logotype

Diff of /MITgcm_contrib/mlosch/interp_llc/addshelfice.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by mlosch, Thu May 3 21:07:20 2007 UTC revision 1.3 by mlosch, Wed May 16 16:09:06 2007 UTC
# Line 1  Line 1 
1  clear all  %clear all
2    
3  new = 'input.180x120x23_shelfice';  new = 'input.shelfice';
4  input = 'input.180x120x23';  input = 'input';
5  eostype = 'mdjwf';  eostype = 'mdjwf';
6    
7  nx = 180;  nx = 45;
8  ny = 120;  ny = nx*18;
9  nz = 23;  nz = 23;
10  nt = 12;  nt = 12;
11    
12    load MASKS
13    hf = msk;
14  load FMT  load FMT
15  hn = mit_readfield(fullfile(input,'bathymetry.bin'),[nx ny],fmt);  load HN
16  hnz = mit_readfield(fullfile(input,'shelfice_bath.bin'),[nx ny],fmt);  load ZN
17  zn = mit_readfield(fullfile(input,'shelfice_topo.bin'),[nx ny],fmt);  % $$$ hn = mit_readfield(fullfile(input,'bathy_llc_p90.bin'),[nx ny],fmt);
18    % $$$ hnz = mit_readfield(fullfile(input,'shelfice_bath.bin'),[nx ny],fmt);
19    % $$$ zn = mit_readfield(fullfile(input,'shelfice_topo.bin'),[nx ny],fmt);
20    
21  h = hn+hnz;  h = hn+hnz;
22  mit_writefield(fullfile(new,'bathymetry.bin'),hn,fmt);  mit_writefield(fullfile(new,'bathy_llc_p90.bin'),mdsiocompact(hn),fmt);
23  mit_writefield(fullfile(new,'bathymetry.shice'),h,fmt);  mit_writefield(fullfile(new,'bathy_llc_p90.shice'),mdsiocompact(h),fmt);
24  mit_writefield(fullfile(new,'shelfice_topo.bin'),zn,fmt);  mit_writefield(fullfile(new,'shelfice_topo.bin'),mdsiocompact(zn),fmt);
25    
26  % create hydrographic fields  % create hydrographic fields
27  levt = mit_readfield(fullfile(input,'lev_t.bin'),[nx ny nz nt],fmt);  levt = mit_readfield(fullfile(input,'lev_t.bin'),[nx ny nz nt],fmt);
# Line 25  levs = mit_readfield(fullfile(input,'lev Line 29  levs = mit_readfield(fullfile(input,'lev
29  is = find(zn~=0);  is = find(zn~=0);
30  [ix,iy] = find(zn~=0);  [ix,iy] = find(zn~=0);
31  [t,s] = shelfice_hydrography(ix,iy,is,levt,levs);  [t,s] = shelfice_hydrography(ix,iy,is,levt,levs);
32  mit_writefield(fullfile(new,'lev_t.shice'),t,fmt);  mit_writefield(fullfile(new,'lev_t.shice'),mdsiocompact(t),fmt);
33  mit_writefield(fullfile(new,'lev_s.shice'),s,fmt);  mit_writefield(fullfile(new,'lev_s.shice'),mdsiocompact(s),fmt);
34    
35    % create hydrographic fields
36    levt = mdsiocompact(mit_readfield(fullfile(input,'lev_t.init'),[nx ny nz],fmt),0);
37    levs = mdsiocompact(mit_readfield(fullfile(input,'lev_s.init'),[nx ny nz],fmt),0);
38    is = find(zn~=0);
39    [ix,iy] = find(zn~=0);
40    [t,s] = shelfice_hydrography(ix,iy,is,levt,levs);
41    mit_writefield(fullfile(new,'lev_t.shice.init'),mdsiocompact(t),fmt);
42    mit_writefield(fullfile(new,'lev_s.shice.init'),mdsiocompact(s),fmt);
43    
44  % create geopotential anomaly  % create geopotential anomaly
45  gravity = 9.81;  gravity = 9.81;
46  rho0 = 1035;  rho0 = 1035;
47  tol = 0;  tol0 = 0;
48  si2dbar = 1e-4;  si2dbar = 1e-4;
 phiHydC = zeros(nz,length(ix));  
 phiHydF = zeros(nz+1,length(ix));  
49  disp('compute geopotential anomaly')  disp('compute geopotential anomaly')
50  load VGRID  load VGRID
51    zg = zf;
52    dzm = abs([zg(1)-zc(1) .5*diff(zc)]);
53    dzp = abs([.5*diff(zc) zc(end)-zg(end)]);
54    hFacMin = 0.1;
55  for ks=1:length(ix)  for ks=1:length(ix)
56    t0 = squeeze(mean(t(ix(ks),iy(ks),:,:),4));    t0 = squeeze(mean(t(ix(ks),iy(ks),:,:),4));
57    s0 = squeeze(mean(s(ix(ks),iy(ks),:,:),4));    s0 = squeeze(mean(s(ix(ks),iy(ks),:,:),4));
58    % compute potential anomaly exactly as in code    % compute potential anomaly exactly as in code
59    % for that we need the correct density    % for that we need the correct density
60    rho = [];    rho = [];
61    p   = -zc(:)*gravity*rho0*si2dbar;    p   = abs(zc(:))*gravity*rho0*si2dbar;
62    dp  = p;    dp  = p;
63    tol1 = 1;    tol1 = 1;
64    tol2 = 2;    tol2 = 2;
65    kw = 0;    kp = 0;
66    while  tol1 > tol    while  tol1 > tol0
67      kw = kw+1;      kp = kp+1;
68        p0 = p;
69      if strcmp(eostype,'mdjwf')      if strcmp(eostype,'mdjwf')
70        rho = [rho densmdjwf(s0,t0,p(:,end))];        drho =  densmdjwf(s0,t0,p(:,end))-rho0;
71      else      else
72        error(['unknown eostype: ' eostype]);        error(['unknown eostype: ' eostype]);
73      end      end
74      p   = [p -zc(:).*rho(:,end)*gravity*si2dbar];      phiHydF(1) = 0;
75        for k=1:length(zc(:));
76          phiHydC(k)   = phiHydF(k) + dzm(k)*gravity*drho(k)/rho0;
77          phiHydF(k+1) = phiHydC(k) + dzp(k)*gravity*drho(k)/rho0;
78        end
79        p   = [p (gravity*rho0*abs(zc(:)) + phiHydC(:)*rho0)/gravity/rho0];
80      dp  = p(:,end)-p(:,end-1);      dp  = p(:,end)-p(:,end-1);
81      tol2 = tol1;      tol2 = tol1;
82      tol1 = sqrt(sum(dp.^2));      tol1 = sqrt(sum(dp.^2));
83      if tol1==tol2; break; end;      if tol1==tol2; break; end;
84    end    end
   % now intergrate  
   drho = rho(:,end)-rho0;  
   for k=1:nz  
     drm = .5*dz(k);  
     if k==1; drm = zf(k)-zc(k); end  
     if k==nz;  
       drp = zc(k)-zf(k+1);  
     else  
       drp = .5*dz(k+1);  
     end  
     phiHydC(k,ks)=phiHydF(k,ks) + drm*gravity*drho(k)/rho0;  
     phiHydF(k+1,ks)=phiHydC(k,ks) +drp*gravity*drho(k)/rho0;  
   end  
85    % find the appropriate level    % find the appropriate level
86    zloc = zn(is(ks));    zloc = zn(is(ks));
87    kl = max(find(zloc < zf));    kl0 = max(find(abs(zg-hFacMin*zg)<=abs(zloc)));
88    ph(ks) = phiHydF(kl,ks);    hfloc= squeeze(hf(ix(ks),iy(ks),:));
89      kl = min(find(hfloc>0));
90      if isempty(kl);
91        kl = 0;
92        ph(ks) = 0;
93      else
94        ph(ks) = phiHydF(kl);
95      end
96      disp(sprintf('kl0 = %u, kl = %u',kl0,kl));
97  end  end
98    
99  pload = zeros(nx,ny);  pload = 0*hn;
100  for ks=1:length(ix)  for ks=1:length(ix)
101    pload(ix(ks),iy(ks)) = -ph(ks)*rho0;    pload(ix(ks),iy(ks)) = -ph(ks)*rho0;
102  end  end
103    
104  mit_writefield(fullfile(new,['pload.' eostype]),pload,fmt);  mit_writefield(fullfile(new,['pload.' eostype]),mdsiocompact(pload),fmt);
105    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22