/[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.2 by mlosch, Fri May 4 14:53:00 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.shelfice';  new = 'input.shelfice';
4  input = 'input';  input = 'input';
# Line 9  ny = nx*18; Line 9  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  load HN  load HN
16  load ZN  load ZN
# Line 42  mit_writefield(fullfile(new,'lev_s.shice Line 44  mit_writefield(fullfile(new,'lev_s.shice
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 = 0*hn;  pload = 0*hn;
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22