1 | function barotropic_currents(clmname,grdname,obc) |
---|
2 | % |
---|
3 | % Pierrick 2003 |
---|
4 | % |
---|
5 | % Get the barotropic velocities from the baroclinic currents |
---|
6 | % Enforce mass conservation |
---|
7 | % |
---|
8 | conserv=1; |
---|
9 | % |
---|
10 | % grid parameters |
---|
11 | % |
---|
12 | disp(' Read grid parameters ...'); |
---|
13 | nc=netcdf(grdname); |
---|
14 | pm=nc{'pm'}(:); |
---|
15 | pn=nc{'pn'}(:); |
---|
16 | h=nc{'h'}(:); |
---|
17 | lon=nc{'lon_rho'}(:); |
---|
18 | lat=nc{'lat_rho'}(:); |
---|
19 | rmask=nc{'mask_rho'}(:); |
---|
20 | umask=nc{'mask_u'}(:); |
---|
21 | vmask=nc{'mask_v'}(:); |
---|
22 | [M,L]=size(rmask); |
---|
23 | close(nc) |
---|
24 | % |
---|
25 | % Model grid vertical levels |
---|
26 | % |
---|
27 | nc=netcdf(clmname,'write'); |
---|
28 | theta_s = nc{'theta_s'}(:); |
---|
29 | theta_b = nc{'theta_b'}(:); |
---|
30 | hc = nc{'hc'}(:); |
---|
31 | N = length(nc('s_rho')); |
---|
32 | tlen = length(nc('uclm_time')); |
---|
33 | % |
---|
34 | % Barotropic velocities |
---|
35 | % |
---|
36 | for l=1:tlen |
---|
37 | disp(['time index: ',num2str(l),' of total: ',num2str(tlen)]) |
---|
38 | zeta=squeeze(nc{'zeta'}(l,:,:)); |
---|
39 | u=squeeze(nc{'u'}(l,:,:,:)); |
---|
40 | v=squeeze(nc{'v'}(l,:,:,:)); |
---|
41 | zw=zlevs(h,zeta,theta_s,theta_b,hc,N,'w'); |
---|
42 | dz=zw(2:end,:,:)-zw(1:end-1,:,:); |
---|
43 | dzu=0.5*(dz(:,:,1:end-1)+dz(:,:,2:end)); |
---|
44 | dzv=0.5*(dz(:,1:end-1,:)+dz(:,2:end,:)); |
---|
45 | hu(:,:)=sum(dzu.*u); |
---|
46 | hv(:,:)=sum(dzv.*v); |
---|
47 | D_u(:,:)=sum(dzu); |
---|
48 | D_v(:,:)=sum(dzv); |
---|
49 | ubar(:,:)=hu./D_u; |
---|
50 | vbar(:,:)=hv./D_v; |
---|
51 | u=u-tridim(ubar,N); |
---|
52 | v=v-tridim(vbar,N); |
---|
53 | % |
---|
54 | % Mass conservation |
---|
55 | % |
---|
56 | if conserv==1 |
---|
57 | disp('Volume conservation enforcement ...') |
---|
58 | [hu,hv]=get_obcvolcons(hu,hv,pm,pn,rmask,obc); |
---|
59 | % |
---|
60 | % Get the stream function |
---|
61 | % |
---|
62 | psi=get_psi(hu,hv,pm,pn,rmask); |
---|
63 | hu(2:end-1,1:end)=-0.5*umask(2:end-1,1:end).*... |
---|
64 | (psi(2:end,1:end)-psi(1:end-1,1:end)).*... |
---|
65 | (pn(2:end-1,2:end)+pn(2:end-1,1:end-1)); |
---|
66 | hv(1:end,2:end-1)=0.5*vmask(1:end,2:end-1).*... |
---|
67 | (psi(1:end,2:end)-psi(1:end,1:end-1)).*... |
---|
68 | (pm(2:end,2:end-1)+pm(1:end-1,2:end-1)); |
---|
69 | [hu,hv]=get_obcvolcons(hu,hv,pm,pn,rmask,obc); |
---|
70 | ubar(:,:)=hu./D_u; |
---|
71 | vbar(:,:)=hv./D_v; |
---|
72 | end |
---|
73 | u=u+tridim(ubar,N); |
---|
74 | v=v+tridim(vbar,N); |
---|
75 | % |
---|
76 | % corners |
---|
77 | % |
---|
78 | ubar(1,1)=0.5*(ubar(1,2)+ubar(2,1)); |
---|
79 | ubar(end,1)=0.5*(ubar(end,2)+ubar(end-1,1)); |
---|
80 | ubar(1,end)=0.5*(ubar(1,end-1)+ubar(2,end)); |
---|
81 | ubar(end,end)=0.5*(ubar(end,end-1)+ubar(end-1,end)); |
---|
82 | vbar(1,1)=0.5*(vbar(1,2)+vbar(2,1)); |
---|
83 | vbar(end,1)=0.5*(vbar(end,2)+vbar(end-1,1)); |
---|
84 | vbar(1,end)=0.5*(vbar(1,end-1)+vbar(2,end)); |
---|
85 | vbar(end,end)=0.5*(vbar(end,end-1)+vbar(end-1,end)); |
---|
86 | u(:,1,1)=0.5*(u(:,1,2)+u(:,2,1)); |
---|
87 | u(:,end,1)=0.5*(u(:,end,2)+u(:,end-1,1)); |
---|
88 | u(:,1,end)=0.5*(u(:,1,end-1)+u(:,2,end)); |
---|
89 | u(:,end,end)=0.5*(u(:,end,end-1)+u(:,end-1,end)); |
---|
90 | v(:,1,1)=0.5*(v(:,1,2)+v(:,2,1)); |
---|
91 | v(:,end,1)=0.5*(v(:,end,2)+v(:,end-1,1)); |
---|
92 | v(:,1,end)=0.5*(v(:,1,end-1)+v(:,2,end)); |
---|
93 | v(:,end,end)=0.5*(v(:,end,end-1)+v(:,end-1,end)); |
---|
94 | % |
---|
95 | % Write into file |
---|
96 | % |
---|
97 | disp('Writes into climatology file ...') |
---|
98 | nc{'u'}(l,:,:,:)=u; |
---|
99 | nc{'v'}(l,:,:,:)=v; |
---|
100 | nc{'ubar'}(l,:,:)=ubar; |
---|
101 | nc{'vbar'}(l,:,:)=vbar; |
---|
102 | end |
---|
103 | close(nc) |
---|