Warning, /verification/solid-body.cs-32x32x1/input/gendata.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit aef96ed3 on 2024-01-05 19:00:20 UTC
aef96ed361 Jean*0001
0002 %-- Generate initial condition (U,V,Eta) for a solid-body rotation flow field
0003 % U-zon = Ueq * cos(Lat) over a sphere in rotation (Omega = 2*pi/rotPeriod).
0004 % Read-in few output grid-files (from dir "gDir") which are easy to get by
0005 % running the MITgcm for 1 time-step with default (=zero) initial conditions.
0006
0007 %- sphere radius [m]:
0008 %rSphere=6370.e+3; %- Earth value
0009 rSphere=5500.4e+3; %- this test
0010 %- rotation period [s] of the solid Earth:
0011 %rotationPeriod= 86400.; % = 1 day
0012 rotationPeriod=108000.; % = 30.h <-- slower rotation rate
0013 %- air density:
0014 rhoConst=1.;
0015 %- Flow velocity (relative to solid Earth) at the Equator [m/s]:
0016 % this relative velocity plus sphere-rotation gives back 24.h rotation period
0017 Ueq=80.;
0018
0019 %- set other constants:
0020 rad=pi/180.;
0021 Omega=2*pi/rotationPeriod; Vsurf=rSphere*Omega;
0022 flowRot=Ueq/rSphere; flowPer= 2*pi/flowRot;
0023
0024 fprintf('Sphere rotation: Omega= %12.6e 1/s ; Period= %6.2f h ; Vsurf= %8.3f m/s\n', ...
0025 Omega, rotationPeriod/3600., Vsurf);
0026 fprintf(' Flow rotation: omega= %12.6e 1/s ; period= %6.2f h ; U.eq = %8.3f m/s\n', ...
0027 flowRot, flowPer/3600., Ueq);
0028
0029 %-- load output grid files fron dir "gDir":
0030 gDir='../run/'; %gDir='./';
0031 xC=rdmds([gDir,'XC']);
0032 yC=rdmds([gDir,'YC']);
0033 yG=rdmds([gDir,'YG']);
0034 dXg=rdmds([gDir,'DXG']);
0035 dYg=rdmds([gDir,'DYG']);
0036
0037 %- put yG into long-vector format and add the 2 missing corners:
0038 n1h=size(dXg,1); n2h=size(dXg,2); nPg=n1h*n2h; nPp2=nPg+2;
0039 yZ=zeros(nPp2,1);
0040 if n1h == 6*n2h,
0041 nc=n2h; %- this the "old-format", all facets stacked along first dimension
0042 yZ(1:nPg,1)=reshape(permute(reshape(yG,[nc 6 nc]),[1 3 2]),[nPg 1]);
0043 yZ(nPg+1)=yG(1+nc*2,1); %- take lat of SW corner of face 3
0044 else
0045 nc=n1h; %- this the "compact-format", all facets stacked along 2nd dimension
0046 yZ(1:nPg,1)=reshape(yG,[nPg 1]);
0047 yZ(nPg+1)=yG(1,1+nc*2); %- take lat of SW corner of face 3
905f660335 Alis*0048 end
aef96ed361 Jean*0049 yZ(nPg+2)=yG(1,1); %- take lat of SW corner of face 1
0050 %- split latitude array in 6 faces:
0051 y6t=split_Z_cub(yZ);
0052
0053 %- define flow field from stream-function:
0054 psi=rSphere*Ueq*sin(rad*y6t);
0055 ut=zeros(nc,nc,6); vt=ut; ncp=nc+1;
0056 for n=1:6,
0057 ut(:,:,n)=psi(1:nc,2:ncp,n)-psi(1:nc,1:nc,n);
0058 vt(:,:,n)=psi(1:nc,1:nc,n)-psi(2:ncp,1:nc,n);
0059 end
0060 if n1h == 6*n2h,
0061 ut=permute(ut,[1 3 2]); vt=permute(vt,[1 3 2]);
0062 end
0063 ut=reshape(ut,[n1h n2h]); vt=reshape(vt,[n1h n2h]);
0064 u=ut./dYg; v=vt./dXg;
0065
0066 %- define surface pressure anomaly from the norm of velocity at grid-cell center
0067 % using analytical expression: cos^2(lat) x fac with:
0068 fac=rhoConst*(Omega*rSphere + Ueq/2)*Ueq;
0069 %- and since shifting by const has no effect on grad(P), substract 2/3(*fac) (= analytical mean)
0070 % to get closer to zero global-mean:
0071 ps1=cos(rad*yC); ps1=ps1.*ps1-2./3.; ps1=fac*ps1;
0072
0073 psAvr=mean(ps1(:));
0074 fprintf('Ps stats: %8.2f , %13.6e , %8.2f\n',min(ps1(:)),psAvr,max(ps1(:)));
0075
0076 %- write initial conditions:
0077 kwr=1; prec='real*8';
0078 % change kwr below to =2 to write (U,V,Eta) initial cond. or to =0 to skip all writing
0079 %kwr=0;
0080 if kwr > 1,
0081 fnam='U_init.bin'; var=u;
0082 fprintf(' writing file: %s ...',fnam);
0083 fid=fopen(fnam,'w','b');
0084 fwrite(fid,var,prec); fclose(fid);
0085 fprintf(' done\n');
0086
0087 fnam='V_init.bin'; var=v;
0088 fprintf(' writing file: %s ...',fnam);
0089 fid=fopen(fnam,'w','b');
0090 fwrite(fid,var,prec); fclose(fid);
0091 fprintf(' done\n');
0092 end
0093 if kwr > 1,
0094 fnam='Eta_ini.bin'; var=ps1;
0095 fprintf(' writing file: %s ...',fnam);
0096 fid=fopen(fnam,'w','b');
0097 fwrite(fid,var,prec); fclose(fid);
0098 fprintf(' done\n');
0099 end
0100
0101 %return
0102
0103 %- define passive tracer (use salt here) patch:
0104 lon=xC*rad; lat=yC*rad;
0105 X=cos(lat).*sin(lon); Y=-cos(lat).*cos(lon); Z=sin(lat);
0106 clear lat lon
0107
0108 lon0=-90 *rad;
0109 lat0=45 *rad;
0110 xo=cos(lat0).*sin(lon0); yo=-cos(lat0).*cos(lon0); zo=sin(lat0);
0111 ro=0.3;
0112 R=sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
0113 s=(1+cos( pi*min(R/ro,1+0*R) ))/2;
0114
0115 if kwr > 0,
0116 fnam='S_init.bin'; var=s;
0117 fprintf(' writing file: %s ...',fnam);
0118 fid=fopen(fnam,'w','b');
0119 fwrite(fid,var,prec); fclose(fid);
0120 fprintf(' done\n');
0121 end
0122
