source: issm/oecreview/Archive/26740-27031/ISSM-26938-26939.diff@ 27032

Last change on this file since 27032 was 27032, checked in by Mathieu Morlighem, 3 years ago

CHG: added 26740-27031
File size: 16.9 KB

  • ../trunk-jpl/test/MITgcm/tools/lookat_4003.m

     
    22p1=[pn 'RunUncoupled/'];
    33p2=[pn 'run/'];
    44
     5for ts=30:34
     6    v1=readbin([p2 'R_shelfIce1_' myint2str(ts,10) '.data'],[3 200]);
     7    v2=readbin([p2 'R_shelfIce2_' myint2str(ts,10) '.data'],[3 200]);
     8    u=readbin([p2 'U.' myint2str(ts,10) '.data'],[3 200 90]);
     9    v=readbin([p2 'V.' myint2str(ts,10) '.data'],[3 200 90]);
     10    w=readbin([p2 'W.' myint2str(ts,10) '.data'],[3 200 90]);
     11    clf, subplot(511), plot(v1(2,:)), title(ts)
     12    subplot(512), plot(v2(2,:)-v1(2,:)), title('draft change')
     13    subplot(513), mypcolor(1:200,-1:-1:-90,squeeze(u(2,:,:))'); title('U'), colorbar
     14    subplot(514), mypcolor(1:200,-1:-1:-90,squeeze(v(2,:,:))'); title('V'), colorbar
     15    subplot(515), mypcolor(1:200,-1:-1:-90,squeeze(w(2,:,:))'); title('W'), colorbar
     16    pause
     17end
     18
     19
     20    clf
     21    subplot(311), mypcolor(v1), title(ts)
     22    subplot(312), mypcolor(v2), title(ts)
     23    subplot(313), mypcolor(v2-v1)
     24   
     25for ts=1:10 
     26    v1=readbin([p2 'SHICE_fwFlux.' myint2str(ts,10) '.data'],[3 200]);
     27    clf
     28    subplot(311), mypcolor(v1), title(ts), colorbar
     29pause
     30end
     31
     32
     33    v2=readbin([p2 'SHICE_fwFlux.' myint2str(ts,10) '.data'],[3 200]);
     34
     35
     36clf
     37for ts=5, disp(ts)
     38    T=readbin([p2 'T.' myint2str(ts,10) '.data'],[3 200 90]);
     39    S=readbin([p2 'S.' myint2str(ts,10) '.data'],[3 200 90]);
     40    U=readbin([p2 'U.' myint2str(ts,10) '.data'],[3 200 90]);
     41    V=readbin([p2 'V.' myint2str(ts,10) '.data'],[3 200 90]);
     42    W=readbin([p2 'W.' myint2str(ts,10) '.data'],[3 200 90]);
     43    for k=1:90, disp(k)
     44        clf
     45        subplot(321), plot(S(2,:,k)), title('S')
     46        subplot(322), plot(T(2,:,k)), title('T')
     47        subplot(323), plot(U(2,:,k)), title('U')
     48        subplot(324), plot(V(2,:,k)), title('V')
     49        subplot(325), plot(W(2,:,k)), title('W')
     50        pause
     51    end
     52end
     53
     54
     55clf
     56for ts=1:8, disp(ts)
     57    T=readbin([p2 'T.' myint2str(ts,10) '.data'],[3 200 90]);
     58    S=readbin([p2 'S.' myint2str(ts,10) '.data'],[3 200 90]);
     59    U=readbin([p2 'U.' myint2str(ts,10) '.data'],[3 200 90]);
     60    V=readbin([p2 'V.' myint2str(ts,10) '.data'],[3 200 90]);
     61    W=readbin([p2 'W.' myint2str(ts,10) '.data'],[3 200 90]);
     62    clf
     63    subplot(321), mypcolor(squeeze(S(2,:,:))'), title('S'), colorbar
     64    subplot(322), mypcolor(squeeze(T(2,:,:))'), title('T'), colorbar
     65    subplot(323), mypcolor(squeeze(U(2,:,:))'), title('U'), colorbar
     66    subplot(324), mypcolor(squeeze(V(2,:,:))'), title('V'), colorbar
     67    subplot(325), mypcolor(squeeze(W(2,:,:))'), title('W'), colorbar
     68    pause
     69end
     70
     71
     72
     73
     74
     75
     76ts=8;
     77v1=readbin([p2 'R_shelfIce1_' myint2str(ts,10) '.data'],[3 200]);
     78v2=readbin([p2 'R_shelfIce2_' myint2str(ts,10) '.data'],[3 200]);
     79clf
     80subplot(311), mypcolor(v1); title('R_shelfIce1'), colorbar
     81subplot(312), mypcolor(v2); title('R_shelfIce2'), colorbar
     82subplot(313), mypcolor(v2-v1); title('diff'), colorbar
     83
     84figure(2)
     85clf
     86plot(1:200,v1(2,:),'o-',1:200,v2(2,:),'o-',1:200,v2(2,:)-v1(2,:),'o-')
     87
    588ts=0;
    689fld='R_shelfIce1_';
    790v1=readbin([p2 fld myint2str(ts,10) '.data'],[3 200]);

  • ../trunk-jpl/test/NightlyRun/test4004.m

     
     1%Test Name: IceOceanCoupling Dan Goldberg'd setup
     2%ISSM/MITgcm coupled set-up
     3%
     4%Script control parameters
     5steps=1:10;
     6
     7%To download and recompile MITgcm from scratch:
     8!rm -rf ${ISSM_DIR}/test/MITgcm/install
     9!rm -rf ${ISSM_DIR}/test/MITgcm/build/*
     10!rm -rf Models
     11
     12%Organizer
     13!mkdir Models
     14org=organizer('repository','Models','prefix','IceOcean.','steps',steps);
     15
     16presentdirectory=pwd;
     17
     18% {{{ Parameters:
     19if perform(org,'Parameters'),
     20    Nx=3;   % number of longitude cells
     21    Ny=200; % number of latitude cells
     22    Nz=90;  % number of MITgcm vertical cells
     23    nPx=1;  % number of MITgcm processes to use in x direction
     24    nPy=8;  % number of MITgcm processes to use in y direction
     25    xgOrigin=0;     % origin of longitude
     26    ygOrigin=-75.5; % origin of latitude
     27    dLong=.125;     % longitude grid spacing
     28    dLat=dLong/16;  % latitude grid spacing
     29    delZ=10;        % thickness of vertical levels (m)
     30    gravity= 9.81;  % gravity (m^2/s)
     31    rho_ice=917;
     32    rho_water=1030;
     33    di=rho_ice/rho_water;
     34    prec = 'real*8'; % precision of MITgcm input binary files
     35
     36    % bathymetry and ice sheet geometry
     37    H = -900;       % water depth in the ice shelf cavity
     38    Hmin = -600;    % deepest point of cavern           
     39    Hmax = -300;    % shallowest point of cavern
     40    jEnd = Ny*3/4;  % where ice-shelf ends
     41    j2   = jEnd+1;
     42   
     43    % initial ocean conditions
     44    T_sfc = -2;
     45    T_bot = -1.9;
     46    S_sfc = 34.2;
     47    S_bot = 34.3;
     48
     49    savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ygOrigin, ...
     50             rho_ice, rho_water, di, H, Hmin, Hmax, jEnd, j2, gravity, prec, ...
     51             T_sfc,T_bot,S_sfc,S_bot);
     52end
     53% }}}
     54% {{{ Bathymetry:
     55if perform(org,'Bathymetry'),
     56
     57    loaddata(org,'Parameters');
     58
     59    %create lat,lon
     60    latg = ygOrigin+[0:Ny-1]*dLat;
     61    latc = latg+dLat/2;
     62    long = xgOrigin+[0:Nx-1]*dLong;
     63    lonc = long+dLong/2;
     64    [lat lon]=meshgrid(latc,lonc);
     65    zC=-delZ*([1:Nz]-0.5);
     66    zF=-delZ*[0:Nz];
     67   
     68    %create bathymetry:
     69    bathymetry = ones(Nx,Ny)*H;
     70    bathymetry(:,end) = 0;
     71
     72    %save bathymetry file for MITgcm
     73    savedata(org,lat,lon,bathymetry);
     74
     75end
     76% }}}
     77% {{{ IceSheetGeometry:
     78if perform(org,'IceSheetGeometry'),
     79
     80    loaddata(org,'Parameters');
     81    loaddata(org,'Bathymetry');
     82    latmin=min(lat(:));
     83    latmax=max(lat(:));
     84
     85    dHdy = (Hmax-Hmin)/dLat/(jEnd-2); %Slope of ice shelf
     86    draft=bathymetry;
     87    for i=1:Nx
     88        draft(i,:)=Hmin+dHdy*[-1:Ny-2]*dLat;
     89    end
     90    draft(:,j2:Ny)=0;
     91
     92    ice_mask=ones(Nx,Ny);
     93    ice_mask(:,j2:Ny)=0;
     94    iceshelf_mask=ice_mask;
     95    thickness=abs(draft)/di;
     96   
     97    savedata(org,ice_mask,iceshelf_mask,draft,thickness);
     98   
     99    close all, figure(2), clf
     100    subplot(411), pcolorcen(bathymetry); colorbar, title('bathymetry')
     101    subplot(412), pcolorcen(ice_mask); colorbar, title('ice and iceshelf mask')
     102    subplot(413), pcolorcen(draft); colorbar, title('draft')
     103    subplot(414), pcolorcen(thickness); colorbar, title('thickness')
     104   
     105end
     106% }}}
     107
     108%Configure MITgcm
     109% {{{ GetMITgcm:
     110if perform(org,'GetMITgcm'),
     111  system([pwd '/../MITgcm/get_mitgcm.sh']);
     112end
     113% }}}
     114% {{{ BuildMITgcm:
     115if perform(org,'BuildMITgcm'),
     116    system(['../MITgcm/build_4003.sh generic ' pwd '/../MITgcm']);
     117end
     118% }}}
     119addpath(recursivepath([pwd '/../MITgcm']));
     120% {{{ RunUncoupledMITgcm:
     121if perform(org,'RunUncoupledMITgcm'),
     122
     123    loaddata(org,'Parameters');
     124    loaddata(org,'Bathymetry');
     125    loaddata(org,'IceSheetGeometry');
     126
     127    % rename previous run directory and create new one
     128    if exist ('run.old')
     129        !\rm -rf run.old
     130    end
     131    if exist ('run')
     132        !\mv run run.old
     133    end
     134    !\mkdir run
     135    !\cp ../MITgcm/build/mitgcmuv run
     136    !\cp ../MITgcm/input_4003/* run
     137    !\cp ../MITgcm/input_4003/eedata_uncoupled run/eedata
     138   
     139    % {{{ Construct MITgcm binary input files
     140    namF='run/bathy_flat.bin';
     141    fid=fopen(namF,'w','b'); fwrite(fid,bathymetry,prec);fclose(fid);
     142   
     143    namF='run/shelficeTopo.Lin.bin';
     144    fid=fopen(namF,'w','b'); fwrite(fid,draft,prec);fclose(fid);
     145   
     146    del_T = (T_bot - T_sfc)/(59*delZ);
     147    tref=zeros(1,Nz);
     148    for k = 1:Nz;
     149        tref(k) = T_sfc + del_T*((k-20)*delZ);
     150        tref(k)= max(T_sfc,min(tref(k),T_bot));
     151    end
     152    namF='run/temp_obc.bin';
     153    tref=[tref; tref; tref];
     154    fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid);
     155
     156    del_S = (S_bot - S_sfc)/(59*delZ);
     157    sref=zeros(1,Nz);
     158    for k = 1:Nz;
     159        sref(k) = S_sfc + del_S*((k-20)*delZ);
     160        sref(k)= max(S_sfc,min(sref(k),S_bot));
     161    end
     162    namF='run/salt_obc.bin';
     163    sref=[sref; sref; sref];
     164    fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid);
     165
     166    zax=[1:Nz];
     167    v1=2.5e-2;
     168    var=1+Nz-2*zax; var=var/(Nz-1);
     169    vobc=v1*var;
     170    namF='run/vVel_obc.bin';
     171    vobc=[vobc; vobc; vobc];
     172    fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid);
     173
     174    var=zeros(Nx,Ny,Nz);
     175    for i=1:Nx, for j=1:Ny
     176            var(i,j,:)=tref(1,:);
     177        end, end
     178    namF='run/temp_ini.bin';
     179    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
     180
     181    for i=1:Nx, for j=1:Ny
     182            var(i,j,:)=sref(1,:);
     183        end, end
     184    namF='run/salt_ini.bin';
     185    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
     186    % }}}
     187
     188    cd run
     189    eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
     190    cd ..
     191end
     192% }}}
     193
     194%Configure ISSM
     195% {{{ CreateMesh:
     196if perform(org,'CreateMesh'),
     197
     198    loaddata(org,'Parameters');
     199
     200    %create model:
     201    md=model();
     202
     203    %Grab lat,long from MITgcm:
     204    long=readbin('run/XC.data',[Nx Ny]);
     205    lat=readbin('run/YC.data',[Nx Ny]);
     206
     207    %project lat,long:
     208    [x,y]=ll2xy(lat(:),long(:),-1);
     209    index=[];
     210    %  C  D
     211    %  A  B
     212    for j=1:(Ny-1),
     213        for i=1:(Nx-1),
     214            A=(j-1)*(Nx)+i;
     215            B=(j-1)*(Nx)+i+1;
     216            C=j*(Nx)+i;
     217            D=j*(Nx)+i+1;
     218            index(end+1,:)=[A B C];
     219            index(end+1,:)=[C B D];
     220        end
     221    end
     222
     223    %fill mesh and model:
     224    md=meshconvert(md,index,x,y);
     225    md.mesh.lat=lat(:);
     226    md.mesh.long=long(:);
     227
     228    savemodel(org,md);
     229
     230end
     231% }}}
     232% {{{ MeshGeometry:
     233if perform(org,'MeshGeometry'),
     234
     235    loaddata(org,'Parameters');
     236    loaddata(org,'CreateMesh');
     237    loaddata(org,'Bathymetry');
     238    loaddata(org,'IceSheetGeometry');
     239
     240    %start filling some of the fields
     241    md.geometry.bed=bathymetry(:);
     242    md.geometry.thickness=thickness(:);
     243    md.geometry.base=-917/1028*md.geometry.thickness;
     244    md.geometry.surface=md.geometry.base+md.geometry.thickness;
     245
     246    %nothing passes icefront:
     247    pos=find((~ice_mask(:) & ice_mask(:)~=0) | thickness(:)==0);
     248    md.geometry.thickness(pos)=1;
     249    md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
     250    md.geometry.base(pos)=-di*md.geometry.thickness(pos);
     251
     252    %level sets:
     253    md.mask.ice_levelset=iceshelf_mask(:);
     254    pos=find(md.mask.ice_levelset==1); md.mask.ice_levelset(pos)=-1;
     255    pos=find(md.mask.ice_levelset==0); md.mask.ice_levelset(pos)=1;
     256    md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
     257
     258    savemodel(org,md);
     259
     260end
     261% }}}
     262% {{{ ParameterizeIce:
     263if perform(org,'ParameterizeIce'),
     264
     265        loaddata(org,'MeshGeometry');
     266
     267        %miscellaneous
     268        md.miscellaneous.name='test4004';
     269
     270        %initial velocity:
     271        md.initialization.vx=zeros(md.mesh.numberofvertices,1);
     272        md.initialization.vy=zeros(md.mesh.numberofvertices,1);
     273        md.initialization.vz=zeros(md.mesh.numberofvertices,1);
     274
     275        %friction:
     276        md.friction.coefficient=0*ones(md.mesh.numberofvertices,1);
     277        pos=find(md.mask.ocean_levelset>0);
     278        md.friction.coefficient(pos)=5;
     279        md.friction.p=ones(md.mesh.numberofelements,1);
     280        md.friction.q=ones(md.mesh.numberofelements,1);
     281
     282        %temperatures and surface mass balance:
     283        md.initialization.temperature=(273.15-22)*ones(md.mesh.numberofvertices,1);
     284        md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
     285        md.smb.mass_balance = 0*ones(md.mesh.numberofvertices,1);
     286
     287        %Flow law
     288        md.materials.rheology_B=paterson(md.initialization.temperature);
     289        md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     290        md.damage.D=zeros(md.mesh.numberofvertices,1);
     291        md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
     292
     293        %the spcs going
     294        md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
     295        md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
     296        md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
     297        md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
     298        md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
     299        md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
     300
     301        %get some flux at the ice divide:
     302        pos=find(md.mesh.lat==min(md.mesh.lat));
     303        md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
     304        md.stressbalance.spcvy(pos)=800;
     305        md.stressbalance.spcvx(pos)=0;
     306
     307        %deal with boundaries, excluding icefront:
     308        pos=find(md.mesh.long==min(md.mesh.long) | md.mesh.long==max(md.mesh.long));
     309        md.stressbalance.spcvx(pos)=0;
     310
     311        point1=find(md.mesh.y==min(md.mesh.y)); point2=find(md.mesh.x==max(md.mesh.x));
     312        costheta=(md.mesh.x(point2)-md.mesh.x(point1))/sqrt((md.mesh.x(point2)-md.mesh.x(point1)).^2+(md.mesh.y(point2)-md.mesh.y(point1)).^2);
     313        sintheta=(md.mesh.y(point2)-md.mesh.y(point1))/sqrt((md.mesh.x(point2)-md.mesh.x(point1)).^2+(md.mesh.y(point2)-md.mesh.y(point1)).^2);
     314        md.stressbalance.referential(:,1:3)=repmat([costheta,sintheta,0],md.mesh.numberofvertices,1);
     315        md.stressbalance.referential(:,4:6)=repmat([-sintheta,costheta,0],md.mesh.numberofvertices,1);
     316
     317        md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     318        md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
     319        md.thermal.spctemperature=[md.initialization.temperature; 1]; %impose observed temperature on surface
     320        md.basalforcings.geothermalflux=.064*ones(md.mesh.numberofvertices,1);
     321
     322        %flow equations:
     323        md=setflowequation(md,'SSA','all');
     324
     325        savemodel(org,md);
     326end
     327% }}}
     328
     329%Run MITgcm/ISSM
     330% {{{ RunCoupledMITgcmISSM:
     331if perform(org,'RunCoupledMITgcmISSM'),
     332
     333    loaddata(org,'Parameters');
     334    loaddata(org,'Bathymetry');
     335    loaddata(org,'IceSheetGeometry');
     336    loaddata(org,'ParameterizeIce');
     337
     338    % {{{ prepare ISSM
     339    md.transient.isoceancoupling=1;
     340    md.transient.isgroundingline=0;
     341    md.transient.isthermal=0;
     342    md.groundingline.migration='SubelementMigration';
     343    md.groundingline.melt_interpolation='NoMeltOnPartiallyFloating';
     344    md.groundingline.friction_interpolation='SubelementFriction2';
     345    md.masstransport.stabilization=1;
     346    md.timestepping.coupling_time=1/24/365/100;          % 1/100 hour in decimal years
     347    md.timestepping.time_step=1/24/365/100;              % 1/100 hour in decimal years
     348    md.timestepping.final_time=23/24/365;            % hour 23 in decimal years
     349    md.cluster=generic('name',oshostname(),'np',2);
     350    md.cluster.npocean=nPx*nPy;
     351    md.cluster.np=2;
     352    md.cluster.executionpath=[pwd '/run'];
     353    md.transient.requested_outputs={'default','MaskOceanLevelset'};
     354    md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};
     355    % }}}
     356    % {{{ prepare MITgcm
     357    % rename previous run directory and create new one
     358    if exist ('run.old')
     359        !\rm -rf run.old
     360    end
     361    if exist ('run')
     362        !\mv run run.old
     363    end
     364    !\mkdir run
     365    !\cp ../MITgcm/build/mitgcmuv run
     366    !\cp ../MITgcm/input_4003/* run
     367   
     368    % {{{ Construct MITgcm binary input files
     369    namF='run/bathy_flat.bin';
     370    fid=fopen(namF,'w','b'); fwrite(fid,bathymetry,prec);fclose(fid);
     371
     372    namF='run/shelficeTopo.Lin.bin';
     373    fid=fopen(namF,'w','b'); fwrite(fid,draft,prec);fclose(fid);
     374
     375    del_T = (T_bot - T_sfc)/(59*delZ);
     376    tref=zeros(1,Nz);
     377    for k = 1:Nz;
     378        tref(k) = T_sfc + del_T*((k-20)*delZ);
     379        tref(k)= max(T_sfc,min(tref(k),T_bot));
     380    end
     381    namF='run/temp_obc.bin';
     382    tref=[tref; tref; tref];
     383    fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid);
     384
     385    del_S = (S_bot - S_sfc)/(59*delZ);
     386    sref=zeros(1,Nz);
     387    for k = 1:Nz;
     388        sref(k) = S_sfc + del_S*((k-20)*delZ);
     389        sref(k)= max(S_sfc,min(sref(k),S_bot));
     390    end
     391    namF='run/salt_obc.bin';
     392    sref=[sref; sref; sref];
     393    fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid);
     394
     395    zax=[1:Nz];
     396    v1=2.5e-2;
     397    var=1+Nz-2*zax; var=var/(Nz-1);
     398    vobc=v1*var;
     399    namF='run/vVel_obc.bin';
     400    vobc=[vobc; vobc; vobc];
     401    fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid);
     402
     403    var=zeros(Nx,Ny,Nz);
     404    for i=1:Nx, for j=1:Ny
     405            var(i,j,:)=tref(1,:);
     406        end, end
     407    namF='run/temp_ini.bin';
     408    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
     409
     410    for i=1:Nx, for j=1:Ny
     411            var(i,j,:)=sref(1,:);
     412        end, end
     413    namF='run/salt_ini.bin';
     414    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
     415    % }}}
     416    % }}}
     417
     418    md=solveiceocean(md,'Transient','runtimename',false);
     419
     420    savemodel(org,md);
     421   
     422    plotmodel(md,'data',md.results.TransientSolution(end).Vel,'data',md.results.TransientSolution(end).Thickness)
     423
     424%eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' ']);
     425%eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' : -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
     426end
     427% }}}
Note: See TracBrowser for help on using the repository browser.