Changeset 24456


Ignore:
Timestamp:
12/07/19 15:54:12 (5 years ago)
Author:
jdquinn
Message:

CHG: Updated config for ice-ocean build; minor cleanup

Location:
issm/trunk-jpl
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/jenkins/ross-debian_linux-iceocean

    r24412 r24456  
    1010        --prefix=$ISSM_DIR \
    1111        --disable-static \
     12        --with-ocean \
    1213        --with-matlab-dir=$MATLAB_PATH \
    1314        --with-fortran-lib="-L/usr/lib/gcc/x86_64-linux-gnu/8 -lgfortran" \
    … …  
    4445        mpich           install-3.3.sh
    4546        lapack          install-3.8-linux.sh
    46         petsc           install-3.11-linux.sh
     47        petsc           install-3.7-linux.sh
    4748        triangle        install-linux64.sh
    4849        chaco           install.sh

  • issm/trunk-jpl/test/NightlyRun/test4001.m

    r24138 r24456  
    66final_time=1;
    77
    8 !rm -rf $ISSM_DIR/test/MITgcm/install
    9 !rm -rf $ISSM_DIR/test/MITgcm/build/*
    10 !rm $ISSM_DIR/test/MITgcm/code/SIZE.h
     8!rm -rf ${ISSM_DIR}/test/MITgcm/install
     9!rm -rf ${ISSM_DIR}/test/MITgcm/build/*
     10!rm -f ${ISSM_DIR}/test/MITgcm/code/SIZE.h
     11!rm -rf Models
    1112
    1213%Organizer
    13 mkdir Models
    14 org=organizer('repository','Models/','prefix','IceOcean.','steps',steps);
    15 
    16 presentdirectory=pwd; 
    17 
    18 % {{{ Parameters: 
     14!mkdir Models
     15org=organizer('repository','Models','prefix','IceOcean.','steps',steps);
     16
     17presentdirectory=pwd;
     18
     19% {{{ Parameters:
    1920if perform(org,'Parameters'),
    20         Nx=20; %number of longitude cells
    21         Ny=40; %number of latitude cells
    22         Nz=30; %number of MITgcm vertical cells
    23         nPx=2; %number of MITgcm processes to use in x direction
    24         nPy=4; %number of MITgcm processes to use in y direction
    25         xgOrigin=0; %origin of longitude
    26         ygOrigin=-80; %origin of latitude
    27         dLong=.25; %longitude grid spacing
    28         dLat=.05; %latitude grid spacing
    29         delZ=30; %thickness of vertical levels
    30         icefront_position_ratio=.75;
    31         ice_thickness=100;
    32         rho_ice=917;
    33         rho_water=1028.14;
    34         di=rho_ice/rho_water;
    35 
    36         % MITgcm initial and lateral boundary conditions
    37         iniSalt  = 34.4; % initial salinity (PSU)
    38         iniTheta = -1.9; % initial potential temperature (deg C)
    39         obcSalt  = 34.4; % open boundary salinity (PSU)
    40         obcTheta =  1.0; % open boundary potential temperature (deg C)
    41         mlDepth  = 120.; % mixed layer depth (m)
    42         mlSalt   = 33.4; % open boundary salinity (PSU)
    43         mlTheta  = -1.9; % open boundary potential temperature (deg C)
    44         obcUvel  = -0.1; % open boundary velocity (m/s)       
    45 
    46         MITgcmDeltaT=600; % MITgcm time step in seconds
    47         y2s=31536000; % year to seconds conversion, i.e., seconds per year
    48 
    49         % start_time, final_time, and time_step
    50         start_time=0; % in decimal years
    51         time_step=1/12; % coupling interval in decimal years
    52         async_step_MITgcm_multiplier=1/30; % used to reduce run time for MItgcm
    53 
    54         % bedrock/bathymetry
    55         hmax=1000;
    56         trough_depth=200;
    57         deltah=300;
    58         sea_level=1095;
    59 
    60         % issm settings:
    61         numlayers=10;
    62 
    63         savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ...
    64                 ygOrigin, icefront_position_ratio, ice_thickness, rho_ice, ...
    65                 rho_water, di, hmax, trough_depth, deltah, sea_level, ...
    66                 iniSalt, iniTheta, obcSalt, obcTheta, mlDepth, mlSalt, ...
    67                 mlTheta, obcUvel, start_time, time_step, MITgcmDeltaT, y2s,...
    68                 numlayers,async_step_MITgcm_multiplier);
    69 end
    70 % }}}
    71 % {{{ Bathymetry: 
     21    Nx=20; %number of longitude cells
     22    Ny=40; %number of latitude cells
     23    Nz=30; %number of MITgcm vertical cells
     24    nPx=2; %number of MITgcm processes to use in x direction
     25    nPy=4; %number of MITgcm processes to use in y direction
     26    xgOrigin=0; %origin of longitude
     27    ygOrigin=-80; %origin of latitude
     28    dLong=.25; %longitude grid spacing
     29    dLat=.05; %latitude grid spacing
     30    delZ=30; %thickness of vertical levels
     31    icefront_position_ratio=.75;
     32    ice_thickness=100;
     33    rho_ice=917;
     34    rho_water=1028.14;
     35    di=rho_ice/rho_water;
     36
     37    % MITgcm initial and lateral boundary conditions
     38    iniSalt  = 34.4; % initial salinity (PSU)
     39    iniTheta = -1.9; % initial potential temperature (deg C)
     40    obcSalt  = 34.4; % open boundary salinity (PSU)
     41    obcTheta =  1.0; % open boundary potential temperature (deg C)
     42    mlDepth  = 120.; % mixed layer depth (m)
     43    mlSalt   = 33.4; % open boundary salinity (PSU)
     44    mlTheta  = -1.9; % open boundary potential temperature (deg C)
     45    obcUvel  = -0.1; % open boundary velocity (m/s)
     46
     47    MITgcmDeltaT=600; % MITgcm time step in seconds
     48    y2s=31536000; % year to seconds conversion, i.e., seconds per year
     49
     50    % start_time, final_time, and time_step
     51    start_time=0; % in decimal years
     52    time_step=1/12; % coupling interval in decimal years
     53    async_step_MITgcm_multiplier=1/30; % used to reduce run time for MItgcm
     54
     55    % bedrock/bathymetry
     56    hmax=1000;
     57    trough_depth=200;
     58    deltah=300;
     59    sea_level=1095;
     60
     61    % issm settings:
     62    numlayers=10;
     63
     64    savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ...
     65        ygOrigin, icefront_position_ratio, ice_thickness, rho_ice, ...
     66        rho_water, di, hmax, trough_depth, deltah, sea_level, ...
     67        iniSalt, iniTheta, obcSalt, obcTheta, mlDepth, mlSalt, ...
     68        mlTheta, obcUvel, start_time, time_step, MITgcmDeltaT, y2s,...
     69        numlayers,async_step_MITgcm_multiplier);
     70end
     71% }}}
     72% {{{ Bathymetry:
    7273if perform(org,'Bathymetry'),
    7374
    … …  
    9495end
    9596% }}}
    96 % {{{ IceSheetGeometry: 
     97% {{{ IceSheetGeometry:
    9798if perform(org,'IceSheetGeometry'),
    98        
    99         loaddata(org,'Parameters');
    100         loaddata(org,'Bathymetry');
    101         latmin=min(lat(:));
    102         latmax=max(lat(:));
    103 
    104         %put ice_thickness constant layer of ice over the bathymetry, unless it floats:
    105         s=size(bathymetry);
    106         thickness=ice_thickness*ones(s);
    107 
    108         %figure out ice shelf:
    109         pos=find(-di*thickness>bathymetry);
    110         iceshelf_mask=zeros(s);
    111         iceshelf_mask(pos)=1;
    112 
    113         ice_mask=ones(s);
    114         pos=find((lat-latmin)/(latmax-latmin)>(icefront_position_ratio));
    115         ice_mask(pos)=0;
    116         iceshelf_mask(pos)=0;
    117 
    118         %compute draft of ice shelf:
    119         draft=bathymetry;
    120         pos=find(iceshelf_mask);
    121         draft(pos)=-di*thickness(pos);
    122         pos=find(~ice_mask);
    123         draft(pos)=0;
    124 
    125         savedata(org,ice_mask,iceshelf_mask,draft,thickness);
     99
     100    loaddata(org,'Parameters');
     101    loaddata(org,'Bathymetry');
     102    latmin=min(lat(:));
     103    latmax=max(lat(:));
     104
     105    %put ice_thickness constant layer of ice over the bathymetry, unless it floats:
     106    s=size(bathymetry);
     107    thickness=ice_thickness*ones(s);
     108
     109    %figure out ice shelf:
     110    pos=find(-di*thickness>bathymetry);
     111    iceshelf_mask=zeros(s);
     112    iceshelf_mask(pos)=1;
     113
     114    ice_mask=ones(s);
     115    pos=find((lat-latmin)/(latmax-latmin)>(icefront_position_ratio));
     116    ice_mask(pos)=0;
     117    iceshelf_mask(pos)=0;
     118
     119    %compute draft of ice shelf:
     120    draft=bathymetry;
     121    pos=find(iceshelf_mask);
     122    draft(pos)=-di*thickness(pos);
     123    pos=find(~ice_mask);
     124    draft(pos)=0;
     125
     126    savedata(org,ice_mask,iceshelf_mask,draft,thickness);
    126127end
    127128% }}}
    128129
    129130%Configure MITgcm
    130 % {{{ GetMITgcm: 
     131% {{{ GetMITgcm:
    131132if perform(org,'GetMITgcm'),
    132133  system([pwd '/../MITgcm/get_mitgcm.sh']);
    133134end
    134135% }}}
    135 % {{{ BuildMITgcm: 
     136% {{{ BuildMITgcm:
    136137if perform(org,'BuildMITgcm'),
    137138
    138     %load data: 
     139    %load data:
    139140    loaddata(org,'Parameters');
    140141
    … …  
    148149            tline = fgetl(fidi);
    149150            if ~ischar(tline), break, end
    150             %do the change here: 
     151            %do the change here:
    151152            if strcmpi(tline,'     &           sNx =  20,'),
    152153                fprintf(fido,'%s%i%s\n','     &           sNx =  ',round(Nx/nPx),',');
    … …  
    176177% }}}
    177178addpath(recursivepath([pwd '/../MITgcm']));
    178 % {{{ RunUncoupledMITgcm: 
     179% {{{ RunUncoupledMITgcm:
    179180if perform(org,'RunUncoupledMITgcm'),
    180181
    181     %load data: 
     182    %load data:
    182183    loaddata(org,'Parameters');
    183184    loaddata(org,'Bathymetry');
    184185    loaddata(org,'IceSheetGeometry');
    185         endtime = round(MITgcmDeltaT * ...
    186                 floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
    187 
    188     % {{{ prepare MITgcm 
     186    endtime = round(MITgcmDeltaT * ...
     187        floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
     188
     189    % {{{ prepare MITgcm
    189190    % rename previous run directory and create new one
    190191    if exist ('run.old')
    … …  
    198199    !\cp ../MITgcm/input/* run
    199200    !\cp ../MITgcm/input/eedata_uncoupled run/eedata
    200    
    201     %load data: 
     201
     202    %load data:
    202203    loaddata(org,'Parameters');
    203204
    … …  
    205206    S=iniSalt*ones(Nx,Ny,Nz);
    206207    writebin('run/Salt.bin',S);
    207    
     208
    208209    % initial temperature
    209210    T=iniTheta*ones(Nx,Ny,Nz);
    210211    writebin('run/Theta.bin',T);
    211    
     212
    212213    % initial velocity
    213214    Z=zeros(Nx,Ny,Nz);
    214215    writebin('run/Uvel.bin',Z);
    215216    writebin('run/Vvel.bin',Z);
    216    
     217
    217218    % initial sea surface height
    218219    Z=zeros(Nx,Ny);
    … …  
    235236    U=obcUvel*ones(Ny,Nz);
    236237    writebin('run/OBu.bin',U);
    237    
     238
    238239    % zero boundary conditions
    239240    Z=zeros(Ny,Nz);
    … …  
    248249        tline = fgetl(fidi);
    249250        if ~ischar(tline), break, end
    250         %do the change here: 
     251        %do the change here:
    251252        if strcmpi(tline,' OB_Iwest = 40*1,'),
    252253            fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
    … …  
    268269    % }}}
    269270
    270     %start looping: 
     271    %start looping:
    271272    for t=start_time:time_step:final_time,
    272273        disp(['Year: ' num2str(t)])
    273         % {{{ generate MITgcm parameter file data 
     274        % {{{ generate MITgcm parameter file data
    274275        fidi=fopen('../MITgcm/input/data','r');
    275276        fido=fopen('run/data','w');
    … …  
    279280            tline = fgetl(fidi);
    280281            if ~ischar(tline), break, end
    281             %do the change here: 
     282            %do the change here:
    282283            if strcmpi(tline,' xgOrigin = 0.0,'),
    283284                fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
    … …  
    369370
    370371%Configure ISSM
    371 % {{{ CreateMesh: 
     372% {{{ CreateMesh:
    372373if perform(org,'CreateMesh'),
    373        
    374         loaddata(org,'Parameters');
    375         loaddata(org,'Bathymetry');
    376         loaddata(org,'IceSheetGeometry');
    377 
    378         %create model:
    379         md=model();
    380        
    381         %Grab lat,long from MITgcm:
    382         lat=lat(:);
    383         long=long(:);
    384 
    385         %project lat,long:
    386         [x,y]=ll2xy(lat,long,-1);
    387 
    388         index=[];
    389         %  C  D
    390         %  A  B
    391         for j=1:Ny-1,
    392                 for i=1:Nx-1,
    393                         A=(j-1)*Nx+i;
    394                         B=(j-1)*Nx+i+1;
    395                         C=j*Nx+i;
    396                         D=j*Nx+i+1;
    397                         index(end+1,:)=[A B C];
    398                         index(end+1,:)=[C B D];
    399                 end
    400         end
    401 
    402         %fill mesh and model:
    403         md=meshconvert(md,index,x,y);
    404         md.mesh.lat=lat;
    405         md.mesh.long=long;
    406 
    407         savemodel(org,md);
    408 
    409 end
    410 % }}}
    411 % {{{ MeshGeometry: 
     374
     375    loaddata(org,'Parameters');
     376    loaddata(org,'Bathymetry');
     377    loaddata(org,'IceSheetGeometry');
     378
     379    %create model:
     380    md=model();
     381
     382    %Grab lat,long from MITgcm:
     383    lat=lat(:);
     384    long=long(:);
     385
     386    %project lat,long:
     387    [x,y]=ll2xy(lat,long,-1);
     388
     389    index=[];
     390    %  C  D
     391    %  A  B
     392    for j=1:Ny-1,
     393        for i=1:Nx-1,
     394            A=(j-1)*Nx+i;
     395            B=(j-1)*Nx+i+1;
     396            C=j*Nx+i;
     397            D=j*Nx+i+1;
     398            index(end+1,:)=[A B C];
     399            index(end+1,:)=[C B D];
     400        end
     401    end
     402
     403    %fill mesh and model:
     404    md=meshconvert(md,index,x,y);
     405    md.mesh.lat=lat;
     406    md.mesh.long=long;
     407
     408    savemodel(org,md);
     409
     410end
     411% }}}
     412% {{{ MeshGeometry:
    412413if perform(org,'MeshGeometry'),
    413        
    414         loaddata(org,'Parameters');
    415         loaddata(org,'CreateMesh');
    416         loaddata(org,'Bathymetry');
    417         loaddata(org,'IceSheetGeometry');
    418 
    419         %transfer to vertices:
    420         bathymetry=bathymetry(:);
    421         iceshelf_mask=iceshelf_mask(:);
    422         ice_mask=ice_mask(:);
    423         thickness=thickness(:);
    424         draft=draft(:);
    425 
    426         %start filling some of the fields
    427         md.geometry.bed=bathymetry;
    428         md.geometry.thickness=thickness;
    429         md.geometry.base=md.geometry.bed;
    430         pos=find(iceshelf_mask); md.geometry.base(pos)=draft(pos);
    431         md.geometry.surface=md.geometry.base+md.geometry.thickness;
    432 
    433         %nothing passes icefront:
    434         pos=find(~ice_mask);
    435         md.geometry.thickness(pos)=1;
    436         md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
    437         md.geometry.base(pos)=-di*md.geometry.thickness(pos);
    438 
    439         %level sets:
    440         md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
    441         md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
    442 
    443         pos=find(ice_mask); md.mask.ice_levelset(pos)=-1;
    444         pos=find(~iceshelf_mask & ice_mask); md.mask.groundedice_levelset(pos)=1;
    445 
    446         %identify edges of grounded ice:
    447         groundedice_levelset=md.mask.groundedice_levelset;
    448         for i=1:md.mesh.numberofelements,
    449                 m=groundedice_levelset(md.mesh.elements(i,:));
    450                 if abs(sum(m))~=3,
    451                         pos=find(m==1); md.mask.groundedice_levelset(md.mesh.elements(i,pos))=0;
    452                 end
    453         end
    454 
    455         %identify edges of ice:
    456         ice_levelset=md.mask.ice_levelset;
    457         for i=1:md.mesh.numberofelements,
    458                 m=ice_levelset(md.mesh.elements(i,:));
    459                 if abs(sum(m))~=3,
    460                         pos=find(m==-1); md.mask.ice_levelset(md.mesh.elements(i,pos))=0;
    461                 end
    462         end
    463 
    464         savemodel(org,md);
    465 end
    466 % }}}
    467 % {{{ ParameterizeIce: 
     414
     415    loaddata(org,'Parameters');
     416    loaddata(org,'CreateMesh');
     417    loaddata(org,'Bathymetry');
     418    loaddata(org,'IceSheetGeometry');
     419
     420    %transfer to vertices:
     421    bathymetry=bathymetry(:);
     422    iceshelf_mask=iceshelf_mask(:);
     423    ice_mask=ice_mask(:);
     424    thickness=thickness(:);
     425    draft=draft(:);
     426
     427    %start filling some of the fields
     428    md.geometry.bed=bathymetry;
     429    md.geometry.thickness=thickness;
     430    md.geometry.base=md.geometry.bed;
     431    pos=find(iceshelf_mask); md.geometry.base(pos)=draft(pos);
     432    md.geometry.surface=md.geometry.base+md.geometry.thickness;
     433
     434    %nothing passes icefront:
     435    pos=find(~ice_mask);
     436    md.geometry.thickness(pos)=1;
     437    md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
     438    md.geometry.base(pos)=-di*md.geometry.thickness(pos);
     439
     440    %level sets:
     441    md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
     442    md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
     443
     444    pos=find(ice_mask); md.mask.ice_levelset(pos)=-1;
     445    pos=find(~iceshelf_mask & ice_mask); md.mask.groundedice_levelset(pos)=1;
     446
     447    %identify edges of grounded ice:
     448    groundedice_levelset=md.mask.groundedice_levelset;
     449    for i=1:md.mesh.numberofelements,
     450        m=groundedice_levelset(md.mesh.elements(i,:));
     451        if abs(sum(m))~=3,
     452            pos=find(m==1); md.mask.groundedice_levelset(md.mesh.elements(i,pos))=0;
     453        end
     454    end
     455
     456    %identify edges of ice:
     457    ice_levelset=md.mask.ice_levelset;
     458    for i=1:md.mesh.numberofelements,
     459        m=ice_levelset(md.mesh.elements(i,:));
     460        if abs(sum(m))~=3,
     461            pos=find(m==-1); md.mask.ice_levelset(md.mesh.elements(i,pos))=0;
     462        end
     463    end
     464
     465    savemodel(org,md);
     466end
     467% }}}
     468% {{{ ParameterizeIce:
    468469if perform(org,'ParameterizeIce'),
    469        
    470         loaddata(org,'Parameters');
    471         loaddata(org,'CreateMesh');
    472         loaddata(org,'MeshGeometry');
    473 
    474         %miscellaneous
    475         md.miscellaneous.name='test4001';
    476 
    477         %initial velocity:
    478         md.initialization.vx=zeros(md.mesh.numberofvertices,1);
    479         md.initialization.vy=zeros(md.mesh.numberofvertices,1);
    480         md.initialization.vz=zeros(md.mesh.numberofvertices,1);
    481 
    482         %friction:
    483         md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
    484         pos=find(md.mask.groundedice_levelset<=0);
    485         md.friction.coefficient(pos)=0;
    486         md.friction.p=ones(md.mesh.numberofelements,1);
    487         md.friction.q=ones(md.mesh.numberofelements,1);
    488 
    489         %temperatures and surface mass balance:
    490         md.initialization.temperature=(273.15-20)*ones(md.mesh.numberofvertices,1);
    491         md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
    492         md.smb.mass_balance = [1*ones(md.mesh.numberofvertices,1); 1];
    493 
    494         %Flow law
    495         md.materials.rheology_B=paterson(md.initialization.temperature);
    496         md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    497         md.damage.D=zeros(md.mesh.numberofvertices,1);
    498         md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
    499        
    500         %the spcs going
    501         md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
    502         md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
    503         md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
    504         md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
    505         md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
    506         md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
    507 
    508         %deal with water:
    509         pos=find(md.mask.ice_levelset>0);
    510         md.stressbalance.spcvx(pos)=0;
    511         md.stressbalance.spcvy(pos)=0;
    512         md.stressbalance.spcvz(pos)=0;
    513         md.masstransport.spcthickness(pos)=0;
    514 
    515         %get some flux at the ice divide:
    516         pos=find(md.mesh.lat==min(md.mesh.lat));
    517         md.stressbalance.spcvy(pos)=200;
    518 
    519         %deal with boundaries, excluding icefront:
    520         vertex_on_boundary=zeros(md.mesh.numberofvertices,1);
    521         vertex_on_boundary(md.mesh.segments(:,1:2))=1;
    522         pos=find(vertex_on_boundary & md.mask.groundedice_levelset<=0);
    523         md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
    524         md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
    525         md.stressbalance.spcvz(pos)=md.initialization.vz(pos);
    526         md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
    527 
    528         md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
    529         md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
    530         md.thermal.spctemperature=[md.initialization.temperature; 1]; %impose observed temperature on surface
    531         md.basalforcings.geothermalflux=.064*ones(md.mesh.numberofvertices,1);
    532 
    533         %flow equations:
    534         md=setflowequation(md,'SSA','all');
    535 
    536         savemodel(org,md);
    537 end
    538 % }}}
    539 % {{{ RunUncoupledISSM: 
     470
     471    loaddata(org,'Parameters');
     472    loaddata(org,'CreateMesh');
     473    loaddata(org,'MeshGeometry');
     474
     475    %miscellaneous
     476    md.miscellaneous.name='test4001';
     477
     478    %initial velocity:
     479    md.initialization.vx=zeros(md.mesh.numberofvertices,1);
     480    md.initialization.vy=zeros(md.mesh.numberofvertices,1);
     481    md.initialization.vz=zeros(md.mesh.numberofvertices,1);
     482
     483    %friction:
     484    md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
     485    pos=find(md.mask.groundedice_levelset<=0);
     486    md.friction.coefficient(pos)=0;
     487    md.friction.p=ones(md.mesh.numberofelements,1);
     488    md.friction.q=ones(md.mesh.numberofelements,1);
     489
     490    %temperatures and surface mass balance:
     491    md.initialization.temperature=(273.15-20)*ones(md.mesh.numberofvertices,1);
     492    md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
     493    md.smb.mass_balance = [1*ones(md.mesh.numberofvertices,1); 1];
     494
     495    %Flow law
     496    md.materials.rheology_B=paterson(md.initialization.temperature);
     497    md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     498    md.damage.D=zeros(md.mesh.numberofvertices,1);
     499    md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
     500
     501    %the spcs going
     502    md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
     503    md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
     504    md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
     505    md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
     506    md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
     507    md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
     508
     509    %deal with water:
     510    pos=find(md.mask.ice_levelset>0);
     511    md.stressbalance.spcvx(pos)=0;
     512    md.stressbalance.spcvy(pos)=0;
     513    md.stressbalance.spcvz(pos)=0;
     514    md.masstransport.spcthickness(pos)=0;
     515
     516    %get some flux at the ice divide:
     517    pos=find(md.mesh.lat==min(md.mesh.lat));
     518    md.stressbalance.spcvy(pos)=200;
     519
     520    %deal with boundaries, excluding icefront:
     521    vertex_on_boundary=zeros(md.mesh.numberofvertices,1);
     522    vertex_on_boundary(md.mesh.segments(:,1:2))=1;
     523    pos=find(vertex_on_boundary & md.mask.groundedice_levelset<=0);
     524    md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
     525    md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
     526    md.stressbalance.spcvz(pos)=md.initialization.vz(pos);
     527    md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
     528
     529    md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     530    md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
     531    md.thermal.spctemperature=[md.initialization.temperature; 1]; %impose observed temperature on surface
     532    md.basalforcings.geothermalflux=.064*ones(md.mesh.numberofvertices,1);
     533
     534    %flow equations:
     535    md=setflowequation(md,'SSA','all');
     536
     537    savemodel(org,md);
     538end
     539% }}}
     540% {{{ RunUncoupledISSM:
    540541if perform(org,'RunUncoupledISSM'),
    541        
    542         loaddata(org,'Parameters');
    543         loaddata(org,'ParameterizeIce');
    544 
    545         %timestepping:
    546         md.timestepping.final_time=final_time;
    547         md.timestepping.time_step=time_step;
    548         md.transient.isgroundingline=1;
    549         md.transient.isthermal=0;
    550         md.groundingline.migration='SubelementMigration';
    551         md.groundingline.melt_interpolation='SubelementMelt2';
    552         md.groundingline.friction_interpolation='SubelementFriction2';
    553 
    554         md.cluster=generic('name',oshostname(),'np',2);
    555         md=solve(md,'Transient');
    556 
    557         savemodel(org,md);
     542
     543    loaddata(org,'Parameters');
     544    loaddata(org,'ParameterizeIce');
     545
     546    %timestepping:
     547    md.timestepping.final_time=final_time;
     548    md.timestepping.time_step=time_step;
     549    md.transient.isgroundingline=1;
     550    md.transient.isthermal=0;
     551    md.groundingline.migration='SubelementMigration';
     552    md.groundingline.melt_interpolation='SubelementMelt2';
     553    md.groundingline.friction_interpolation='SubelementFriction2';
     554
     555    md.cluster=generic('name',oshostname(),'np',2);
     556    md=solve(md,'Transient');
     557
     558    savemodel(org,md);
    558559end
    559560% }}}
    560561
    561562%Run MITgcm/ISSM
    562 % {{{ RunCoupledMITgcmISSM: 
     563% {{{ RunCoupledMITgcmISSM:
    563564if perform(org,'RunCoupledMITgcmISSM'),
    564565
    565         %load data:
    566         loaddata(org,'Parameters');
    567         loaddata(org,'ParameterizeIce');
    568         loaddata(org,'Bathymetry');
    569         loaddata(org,'IceSheetGeometry');
     566    %load data:
     567    loaddata(org,'Parameters');
     568    loaddata(org,'ParameterizeIce');
     569    loaddata(org,'Bathymetry');
     570    loaddata(org,'IceSheetGeometry');
    570571        endtime = round(MITgcmDeltaT * ...
    571572         floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
    572573
    573         % {{{ prepare MITgcm 
     574        % {{{ prepare MITgcm
    574575        % rename previous run directory and create new one
    575576        if exist ('run.old')
    … …  
    584585        !\cp ../MITgcm/input/eedata_uncoupled run/eedata
    585586
    586         %load data: 
     587        %load data:
    587588        loaddata(org,'Parameters');
    588589
    … …  
    633634            tline = fgetl(fidi);
    634635            if ~ischar(tline), break, end
    635             %do the change here: 
     636            %do the change here:
    636637            if strcmpi(tline,' OB_Iwest = 40*1,'),
    637638                fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
    … …  
    651652        writebin('run/bathymetry.bin',bathymetry);
    652653        % }}}
    653            
    654         % {{{ ISSM settings:
    655 
    656         setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib')
    657         %timestepping:
    658         md.timestepping.start_time=start_time;
    659         md.timestepping.final_time=final_time;
    660         md.timestepping.time_step=time_step;
    661         md.cluster=generic('name',oshostname(),'np',2);
    662         md.results.TransientSolution.Base=md.geometry.base;
    663         md.transient.isgroundingline=1;
    664         md.transient.isthermal=0;
    665         md.groundingline.migration='SubelementMigration';
    666         md.groundingline.melt_interpolation='SubelementMelt2';
    667         md.groundingline.friction_interpolation='SubelementFriction2';
    668 
    669         % }}}
    670 
    671         %start looping:
    672         results=md.results;
    673 
    674         for t=start_time:time_step:final_time
     654
     655    % {{{ ISSM settings:
     656
     657    setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib')
     658    %timestepping:
     659    md.timestepping.start_time=start_time;
     660    md.timestepping.final_time=final_time;
     661    md.timestepping.time_step=time_step;
     662    md.cluster=generic('name',oshostname(),'np',2);
     663    md.results.TransientSolution.Base=md.geometry.base;
     664    md.transient.isgroundingline=1;
     665    md.transient.isthermal=0;
     666    md.groundingline.migration='SubelementMigration';
     667    md.groundingline.melt_interpolation='SubelementMelt2';
     668    md.groundingline.friction_interpolation='SubelementFriction2';
     669
     670    % }}}
     671
     672    %start looping:
     673    results=md.results;
     674
     675    for t=start_time:time_step:final_time
    675676            disp(['Year: ' num2str(t)])
    676677
    677             %send draft from ISSM to MITgcm:
    678             draft=md.results.TransientSolution(end).Base;
    679             pos=find(md.mask.ice_levelset>0); draft(pos)=0;
     678        %send draft from ISSM to MITgcm:
     679        draft=md.results.TransientSolution(end).Base;
     680        pos=find(md.mask.ice_levelset>0); draft(pos)=0;
    680681            if t>start_time
    681682                old_draft=readbin('run/icetopo.bin',[Nx,Ny]);
    … …  
    683684            writebin('run/icetopo.bin',draft);
    684685
    685             % {{{ generate MITgcm parameter file data
    686             fidi=fopen('../MITgcm/input/data','r');
    687             fido=fopen('run/data','w');
    688             tline = fgetl(fidi);
    689             fprintf(fido,'%s\n',tline);
     686        % {{{ generate MITgcm parameter file data
     687        fidi=fopen('../MITgcm/input/data','r');
     688        fido=fopen('run/data','w');
     689        tline = fgetl(fidi);
     690        fprintf(fido,'%s\n',tline);
    690691            while 1
    691692                tline = fgetl(fidi);
    692693                if ~ischar(tline), break, end
    693                 %do the change here:
    694                 if strcmpi(tline,' xgOrigin = 0.0,'),
    695                     fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
    696                     continue;
    697                 end
    698                 if strcmpi(tline,' ygOrigin = -80.0,'),
    699                     fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
    700                     continue;
    701                 end
    702                 if strcmpi(tline,' delX = 20*0.25,'),
    703                     fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
    704                     continue;
    705                 end
    706                 if strcmpi(tline,' delY = 20*0.25,'),
    707                     fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
    708                     continue;
    709                 end
     694        %do the change here:
     695        if strcmpi(tline,' xgOrigin = 0.0,'),
     696            fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
     697            continue;
     698        end
     699        if strcmpi(tline,' ygOrigin = -80.0,'),
     700            fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
     701            continue;
     702        end
     703        if strcmpi(tline,' delX = 20*0.25,'),
     704            fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
     705            continue;
     706        end
     707        if strcmpi(tline,' delY = 20*0.25,'),
     708            fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
     709            continue;
     710        end
    710711                if strcmpi(tline,' delZ = 30*30.0,'),
    711712                    fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
    712713                    continue;
    713714                end
    714                 if strcmpi(tline,' endTime=2592000.,'),
    715                     fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
    716                     continue;
    717                 end
    718                 if strcmpi(tline,' deltaT=1200.0,'),
    719                     fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
    720                     continue;
    721                 end
    722                 if strcmpi(tline,' pChkptFreq=2592000.,'),
    723                     fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
    724                     continue;
    725                 end
    726                 if strcmpi(tline,' taveFreq=2592000.,'),
    727                     fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
    728                     continue;
    729                 end
     715        if strcmpi(tline,' endTime=2592000.,'),
     716            fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
     717            continue;
     718        end
     719        if strcmpi(tline,' deltaT=1200.0,'),
     720            fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
     721            continue;
     722        end
     723        if strcmpi(tline,' pChkptFreq=2592000.,'),
     724            fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
     725            continue;
     726        end
     727        if strcmpi(tline,' taveFreq=2592000.,'),
     728            fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
     729            continue;
     730        end
    730731                if strcmpi(tline,' rhoConst=1030.,'),
    731732                    fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
    … …  
    736737                    continue;
    737738                end
    738                 fprintf(fido,'%s\n',tline);
    739             end
    740             %close  files
    741             fclose(fidi);
    742             fclose(fido);
    743             % }}}
    744 
    745             % {{{ generate initial MITgcm conditions
     739        fprintf(fido,'%s\n',tline);
     740        end
     741        %close  files
     742        fclose(fidi);
     743        fclose(fido);
     744        % }}}
     745
     746        % {{{ generate initial MITgcm conditions
    746747            ds=round(endtime/MITgcmDeltaT);
    747748            if t>start_time
    … …  
    762763                tmp(find(tmp<0))=0;
    763764                [im jm]=find(tmp); % horizontal indices where there is melt
    764                
     765
    765766                % Extrapolate T/S to locations where ice shelf retreated
    766767                for i=1:length(im)
    … …  
    796797            % }}}
    797798
    798             % {{{ system call to run MITgcm 
     799            % {{{ system call to run MITgcm
    799800            cd run
    800801            eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
    … …  
    812813            % }}}
    813814
    814             %get melting rates from MITgcm
    815             %upward fresh water flux (kg/m^2/s):
    816             fnm=['run/SHICE_fwFluxtave_' myint2str(ts,10) '.data'];
    817             melting_rate=readbin(fnm,[Nx Ny]);
    818 
    819             %send averaged melting rate to ISSM
    820             %downward fresh water flux (m/y):
    821             melting_rate=-melting_rate(:)*y2s/rho_ice;
    822                 md.basalforcings.floatingice_melting_rate=melting_rate;
    823 
    824             % {{{ run ISSM and recover results
    825 
    826             md.timestepping.start_time=t;
    827             md.timestepping.final_time=t+time_step;;
    828                 md=solve(md,'Transient');
    829 
    830                 base=md.results.TransientSolution(end).Base;
    831                 thickness=md.results.TransientSolution(end).Thickness;
    832                 md.geometry.base=base;
    833                 md.geometry.thickness=thickness;
    834                 md.geometry.surface=md.geometry.base+md.geometry.thickness;
    835                 md.initialization.vx=md.results.TransientSolution(end).Vx;
    836                 md.initialization.vy=md.results.TransientSolution(end).Vy;
    837                 md.initialization.vel=md.results.TransientSolution(end).Vel;
    838                 md.initialization.pressure=md.results.TransientSolution(end).Pressure;
    839                 md.mask.groundedice_levelset=md.results.TransientSolution(end).MaskGroundediceLevelset;
    840                 md.results.TransientSolution(end).FloatingiceMeltingRate=md.basalforcings.floatingice_melting_rate;
    841                
    842                 %save these results in the model, otherwise, they'll be wiped out
    843                 results(end+1)=md.results;
    844 
    845                 % }}}
    846 
    847         end
    848 
    849         md.results=results;
    850         savemodel(org,md);
     815        %get melting rates from MITgcm
     816        %upward fresh water flux (kg/m^2/s):
     817        fnm=['run/SHICE_fwFluxtave_' myint2str(ts,10) '.data'];
     818        melting_rate=readbin(fnm,[Nx Ny]);
     819
     820        %send averaged melting rate to ISSM
     821        %downward fresh water flux (m/y):
     822        melting_rate=-melting_rate(:)*y2s/rho_ice;
     823        md.basalforcings.floatingice_melting_rate=melting_rate;
     824
     825        % {{{ run ISSM and recover results
     826
     827        md.timestepping.start_time=t;
     828        md.timestepping.final_time=t+time_step;;
     829        md=solve(md,'Transient');
     830
     831        base=md.results.TransientSolution(end).Base;
     832        thickness=md.results.TransientSolution(end).Thickness;
     833        md.geometry.base=base;
     834        md.geometry.thickness=thickness;
     835        md.geometry.surface=md.geometry.base+md.geometry.thickness;
     836        md.initialization.vx=md.results.TransientSolution(end).Vx;
     837        md.initialization.vy=md.results.TransientSolution(end).Vy;
     838        md.initialization.vel=md.results.TransientSolution(end).Vel;
     839        md.initialization.pressure=md.results.TransientSolution(end).Pressure;
     840        md.mask.groundedice_levelset=md.results.TransientSolution(end).MaskGroundediceLevelset;
     841        md.results.TransientSolution(end).FloatingiceMeltingRate=md.basalforcings.floatingice_melting_rate;
     842
     843        %save these results in the model, otherwise, they'll be wiped out
     844        results(end+1)=md.results;
     845
     846        % }}}
     847
     848    end
     849
     850    md.results=results;
     851    savemodel(org,md);
    851852end
    852853% }}}
    … …  
    862863melting_rate_4=readbin(fnm,[Nx Ny]);
    863864field_names     ={'Base1','Melting1','Vx2','Vy2','Thickness2','Base2','MaskGroundediceLevelset2','FloatingiceMeltingRate2',...
    864         'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',...
    865         'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'};
     865    'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',...
     866    'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'};
    866867field_tolerances={2e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,...
    867         1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,...
    868         1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 };
     868    1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,...
     869    1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 };
    869870field_values={...
    870         (md.results(1).TransientSolution(end).Base),...
    871         (melting_rate_1(:)),...
    872         (md.results(2).TransientSolution(end).Vx),...
    873         (md.results(2).TransientSolution(end).Vy),...
    874         (md.results(2).TransientSolution(end).Thickness),...
    875         (md.results(2).TransientSolution(end).Base),...
    876         (md.results(2).TransientSolution(end).MaskGroundediceLevelset),...
    877         (md.results(2).TransientSolution(end).FloatingiceMeltingRate),...
    878         (melting_rate_2(:)),...
    879         (md.results(3).TransientSolution(end).Vx),...
    880         (md.results(3).TransientSolution(end).Vy),...
    881         (md.results(3).TransientSolution(end).Thickness),...
    882         (md.results(3).TransientSolution(end).Base),...
    883         (md.results(3).TransientSolution(end).MaskGroundediceLevelset),...
    884         (md.results(3).TransientSolution(end).FloatingiceMeltingRate),...
    885         (melting_rate_3(:)),...
    886         (md.results(4).TransientSolution(end).Vx),...
    887         (md.results(4).TransientSolution(end).Vy),...
    888         (md.results(4).TransientSolution(end).Thickness),...
    889         (md.results(4).TransientSolution(end).Base),...
    890         (md.results(4).TransientSolution(end).MaskGroundediceLevelset),...
    891         (md.results(4).TransientSolution(end).FloatingiceMeltingRate),...
    892         (melting_rate_4(:)),...
    893         };
     871    (md.results(1).TransientSolution(end).Base),...
     872    (melting_rate_1(:)),...
     873    (md.results(2).TransientSolution(end).Vx),...
     874    (md.results(2).TransientSolution(end).Vy),...
     875    (md.results(2).TransientSolution(end).Thickness),...
     876    (md.results(2).TransientSolution(end).Base),...
     877    (md.results(2).TransientSolution(end).MaskGroundediceLevelset),...
     878    (md.results(2).TransientSolution(end).FloatingiceMeltingRate),...
     879    (melting_rate_2(:)),...
     880    (md.results(3).TransientSolution(end).Vx),...
     881    (md.results(3).TransientSolution(end).Vy),...
     882    (md.results(3).TransientSolution(end).Thickness),...
     883    (md.results(3).TransientSolution(end).Base),...
     884    (md.results(3).TransientSolution(end).MaskGroundediceLevelset),...
     885    (md.results(3).TransientSolution(end).FloatingiceMeltingRate),...
     886    (melting_rate_3(:)),...
     887    (md.results(4).TransientSolution(end).Vx),...
     888    (md.results(4).TransientSolution(end).Vy),...
     889    (md.results(4).TransientSolution(end).Thickness),...
     890    (md.results(4).TransientSolution(end).Base),...
     891    (md.results(4).TransientSolution(end).MaskGroundediceLevelset),...
     892    (md.results(4).TransientSolution(end).FloatingiceMeltingRate),...
     893    (melting_rate_4(:)),...
     894    };

  • issm/trunk-jpl/test/NightlyRun/test4002.m

    r23881 r24456  
    33%
    44%Script control parameters
    5 steps=[1 2 3 4 5 6 7 8 9 10 11 12];
     5steps=1:12;
    66final_time=1/365;
    77
    88%To download and recompile MITgcm from scratch:
    9 !rm -rf $ISSM_DIR/test/MITgcm/install
    10 !rm -rf $ISSM_DIR/test/MITgcm/build/*
    11 !rm $ISSM_DIR/test/MITgcm/code/SIZE.h
     9!rm -rf ${ISSM_DIR}/test/MITgcm/install
     10!rm -rf ${ISSM_DIR}/test/MITgcm/build/*
     11!rm -f ${ISSM_DIR}/test/MITgcm/code/SIZE.h
     12!rm -rf Models
    1213
    1314%Organizer
    14 mkdir Models
    15 org=organizer('repository','Models/','prefix','IceOcean.','steps',steps);
     15!mkdir Models
     16org=organizer('repository','Models','prefix','IceOcean.','steps',steps);
    1617
    1718presentdirectory=pwd;
    1819
    19 % {{{ Parameters: 
     20% {{{ Parameters:
    2021if perform(org,'Parameters'),
    21         Nx=20; %number of longitude cells
    22         Ny=40; %number of latitude cells
    23         Nz=30; %number of MITgcm vertical cells
    24         nPx=2; %number of MITgcm processes to use in x direction
    25         nPy=4; %number of MITgcm processes to use in y direction
    26         xgOrigin=0; %origin of longitude
    27         ygOrigin=-80; %origin of latitude
    28         dLong=.25; %longitude grid spacing
    29         dLat=.05; %latitude grid spacing
    30         delZ=30; %thickness of vertical levels
    31         icefront_position_ratio=.75;
    32         ice_thickness=100;
    33         rho_ice=917;
    34         rho_water=1028.14;
    35         di=rho_ice/rho_water;
    36 
    37         % MITgcm initial and lateral boundary conditions
    38         iniSalt  = 34.4; % initial salinity (PSU)
    39         iniTheta = -1.9; % initial potential temperature (deg C)
    40         obcSalt  = 34.4; % open boundary salinity (PSU)
    41         obcTheta =  1.0; % open boundary potential temperature (deg C)
    42         mlDepth  = 120.; % mixed layer depth (m)
    43         mlSalt   = 33.4; % open boundary salinity (PSU)
    44         mlTheta  = -1.9; % open boundary potential temperature (deg C)
    45         obcUvel  = -0.1; % open boundary velocity (m/s)       
    46 
    47         MITgcmDeltaT=600; % MITgcm time step in seconds
    48         y2s=31536000; % year to seconds conversion, i.e., seconds per year
    49 
    50         % start_time and time_step
    51         start_time=0; % in decimal years
    52         time_step=1/(365*24); % coupling interval in decimal years
    53         async_step_MITgcm_multiplier=1; % used to reduce run time for MItgcm
    54 
    55         % bedrock/bathymetry
    56         hmax=1000;
    57         trough_depth=200;
    58         deltah=300;
    59         sea_level=1095;
    60 
    61         % issm settings:
    62         numlayers=10;
    63 
    64         savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ...
    65                 ygOrigin, icefront_position_ratio, ice_thickness, rho_ice, ...
    66                 rho_water, di, hmax, trough_depth, deltah, sea_level, ...
    67                 iniSalt, iniTheta, obcSalt, obcTheta, mlDepth, mlSalt, ...
    68                 mlTheta, obcUvel, start_time, time_step, MITgcmDeltaT, y2s,...
    69                 numlayers,async_step_MITgcm_multiplier);
    70 end
    71 % }}}
    72 % {{{ Bathymetry: 
     22    Nx=20; %number of longitude cells
     23    Ny=40; %number of latitude cells
     24    Nz=30; %number of MITgcm vertical cells
     25    nPx=2; %number of MITgcm processes to use in x direction
     26    nPy=4; %number of MITgcm processes to use in y direction
     27    xgOrigin=0; %origin of longitude
     28    ygOrigin=-80; %origin of latitude
     29    dLong=.25; %longitude grid spacing
     30    dLat=.05; %latitude grid spacing
     31    delZ=30; %thickness of vertical levels
     32    icefront_position_ratio=.75;
     33    ice_thickness=100;
     34    rho_ice=917;
     35    rho_water=1028.14;
     36    di=rho_ice/rho_water;
     37
     38    % MITgcm initial and lateral boundary conditions
     39    iniSalt  = 34.4; % initial salinity (PSU)
     40    iniTheta = -1.9; % initial potential temperature (deg C)
     41    obcSalt  = 34.4; % open boundary salinity (PSU)
     42    obcTheta =  1.0; % open boundary potential temperature (deg C)
     43    mlDepth  = 120.; % mixed layer depth (m)
     44    mlSalt   = 33.4; % open boundary salinity (PSU)
     45    mlTheta  = -1.9; % open boundary potential temperature (deg C)
     46    obcUvel  = -0.1; % open boundary velocity (m/s)
     47
     48    MITgcmDeltaT=600; % MITgcm time step in seconds
     49    y2s=31536000; % year to seconds conversion, i.e., seconds per year
     50
     51    % start_time and time_step
     52    start_time=0; % in decimal years
     53    time_step=1/(365*24); % coupling interval in decimal years
     54    async_step_MITgcm_multiplier=1; % used to reduce run time for MItgcm
     55
     56    % bedrock/bathymetry
     57    hmax=1000;
     58    trough_depth=200;
     59    deltah=300;
     60    sea_level=1095;
     61
     62    % issm settings:
     63    numlayers=10;
     64
     65    savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ...
     66        ygOrigin, icefront_position_ratio, ice_thickness, rho_ice, ...
     67        rho_water, di, hmax, trough_depth, deltah, sea_level, ...
     68        iniSalt, iniTheta, obcSalt, obcTheta, mlDepth, mlSalt, ...
     69        mlTheta, obcUvel, start_time, time_step, MITgcmDeltaT, y2s,...
     70        numlayers,async_step_MITgcm_multiplier);
     71end
     72% }}}
     73% {{{ Bathymetry:
    7374if perform(org,'Bathymetry'),
    7475
    … …  
    9596end
    9697% }}}
    97 % {{{ IceSheetGeometry: 
     98% {{{ IceSheetGeometry:
    9899if perform(org,'IceSheetGeometry'),
    99        
    100         loaddata(org,'Parameters');
    101         loaddata(org,'Bathymetry');
    102         latmin=min(lat(:));
    103         latmax=max(lat(:));
    104 
    105         %put ice_thickness constant layer of ice over the bathymetry, unless it floats:
    106         s=size(bathymetry);
    107         thickness=ice_thickness*ones(s);
    108 
    109         %figure out ice shelf:
    110         pos=find(-di*thickness>bathymetry);
    111         iceshelf_mask=zeros(s);
    112         iceshelf_mask(pos)=1;
    113 
    114         ice_mask=ones(s);
    115         pos=find((lat-latmin)/(latmax-latmin)>(icefront_position_ratio));
    116         ice_mask(pos)=0;
    117         iceshelf_mask(pos)=0;
    118 
    119         %compute draft of ice shelf:
    120         draft=bathymetry;
    121         pos=find(iceshelf_mask);
    122         draft(pos)=-di*thickness(pos);
    123         pos=find(~ice_mask);
    124         draft(pos)=0;
    125 
    126         savedata(org,ice_mask,iceshelf_mask,draft,thickness);
     100
     101    loaddata(org,'Parameters');
     102    loaddata(org,'Bathymetry');
     103    latmin=min(lat(:));
     104    latmax=max(lat(:));
     105
     106    %put ice_thickness constant layer of ice over the bathymetry, unless it floats:
     107    s=size(bathymetry);
     108    thickness=ice_thickness*ones(s);
     109
     110    %figure out ice shelf:
     111    pos=find(-di*thickness>bathymetry);
     112    iceshelf_mask=zeros(s);
     113    iceshelf_mask(pos)=1;
     114
     115    ice_mask=ones(s);
     116    pos=find((lat-latmin)/(latmax-latmin)>(icefront_position_ratio));
     117    ice_mask(pos)=0;
     118    iceshelf_mask(pos)=0;
     119
     120    %compute draft of ice shelf:
     121    draft=bathymetry;
     122    pos=find(iceshelf_mask);
     123    draft(pos)=-di*thickness(pos);
     124    pos=find(~ice_mask);
     125    draft(pos)=0;
     126
     127    savedata(org,ice_mask,iceshelf_mask,draft,thickness);
    127128end
    128129% }}}
    129130
    130131%Configure MITgcm
    131 % {{{ GetMITgcm: 
     132% {{{ GetMITgcm:
    132133if perform(org,'GetMITgcm'),
    133134  system([pwd '/../MITgcm/get_mitgcm.sh']);
    134135end
    135136% }}}
    136 % {{{ BuildMITgcm: 
     137% {{{ BuildMITgcm:
    137138if perform(org,'BuildMITgcm'),
    138139
    139     %load data: 
     140    %load data:
    140141    loaddata(org,'Parameters');
    141142
    … …  
    149150            tline = fgetl(fidi);
    150151            if ~ischar(tline), break, end
    151             %do the change here: 
     152            %do the change here:
    152153            if strcmpi(tline,'     &           sNx =  20,'),
    153154                fprintf(fido,'%s%i%s\n','     &           sNx =  ',round(Nx/nPx),',');
    … …  
    177178% }}}
    178179addpath(recursivepath([pwd '/../MITgcm']));
    179 % {{{ RunUncoupledMITgcm: 
     180% {{{ RunUncoupledMITgcm:
    180181if perform(org,'RunUncoupledMITgcm'),
    181182
    182     %load data: 
     183    %load data:
    183184    loaddata(org,'Parameters');
    184185    loaddata(org,'Bathymetry');
    185186    loaddata(org,'IceSheetGeometry');
    186         endtime = round(MITgcmDeltaT * ...
    187                 floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
    188 
    189     % {{{ prepare MITgcm 
     187    endtime = round(MITgcmDeltaT * ...
     188        floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
     189
     190    % {{{ prepare MITgcm
    190191    % rename previous run directory and create new one
    191192    if exist ('run.old')
    … …  
    199200    !\cp ../MITgcm/input/* run
    200201    !\cp ../MITgcm/input/eedata_uncoupled run/eedata
    201    
    202     %load data: 
     202
     203    %load data:
    203204    loaddata(org,'Parameters');
    204205
    … …  
    206207    S=iniSalt*ones(Nx,Ny,Nz);
    207208    writebin('run/Salt.bin',S);
    208    
     209
    209210    % initial temperature
    210211    T=iniTheta*ones(Nx,Ny,Nz);
    211212    writebin('run/Theta.bin',T);
    212    
     213
    213214    % initial velocity
    214215    Z=zeros(Nx,Ny,Nz);
    215216    writebin('run/Uvel.bin',Z);
    216217    writebin('run/Vvel.bin',Z);
    217    
     218
    218219    % initial sea surface height
    219220    Z=zeros(Nx,Ny);
    … …  
    236237    U=obcUvel*ones(Ny,Nz);
    237238    writebin('run/OBu.bin',U);
    238    
     239
    239240    % zero boundary conditions
    240241    Z=zeros(Ny,Nz);
    … …  
    249250        tline = fgetl(fidi);
    250251        if ~ischar(tline), break, end
    251         %do the change here: 
     252        %do the change here:
    252253        if strcmpi(tline,' OB_Iwest = 40*1,'),
    253254            fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
    … …  
    269270    % }}}
    270271
    271     %start looping: 
     272    %start looping:
    272273    for t=start_time:time_step:final_time,
    273274        disp(['Year: ' num2str(t)])
    274         % {{{ generate MITgcm parameter file data 
     275        % {{{ generate MITgcm parameter file data
    275276        fidi=fopen('../MITgcm/input/data','r');
    276277        fido=fopen('run/data','w');
    … …  
    280281            tline = fgetl(fidi);
    281282            if ~ischar(tline), break, end
    282             %do the change here: 
     283            %do the change here:
    283284            if strcmpi(tline,' xgOrigin = 0.0,'),
    284285                fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
    … …  
    370371
    371372%Configure ISSM
    372 % {{{ CreateMesh: 
     373% {{{ CreateMesh:
    373374if perform(org,'CreateMesh'),
    374        
    375         loaddata(org,'Parameters');
    376         loaddata(org,'Bathymetry');
    377         loaddata(org,'IceSheetGeometry');
    378 
    379         %create model:
    380         md=model();
    381        
    382         %Grab lat,long from MITgcm:
    383         lat=lat(:);
    384         long=long(:);
    385 
    386         %project lat,long:
    387         [x,y]=ll2xy(lat,long,-1);
    388 
    389         index=[];
    390         %  C  D
    391         %  A  B
    392         for j=1:Ny-1,
    393                 for i=1:Nx-1,
    394                         A=(j-1)*Nx+i;
    395                         B=(j-1)*Nx+i+1;
    396                         C=j*Nx+i;
    397                         D=j*Nx+i+1;
    398                         index(end+1,:)=[A B C];
    399                         index(end+1,:)=[C B D];
    400                 end
    401         end
    402 
    403         %fill mesh and model:
    404         md=meshconvert(md,index,x,y);
    405         md.mesh.lat=lat;
    406         md.mesh.long=long;
    407 
    408         savemodel(org,md);
    409 
    410 end
    411 % }}}
    412 % {{{ MeshGeometry: 
     375
     376    loaddata(org,'Parameters');
     377    loaddata(org,'Bathymetry');
     378    loaddata(org,'IceSheetGeometry');
     379
     380    %create model:
     381    md=model();
     382
     383    %Grab lat,long from MITgcm:
     384    lat=lat(:);
     385    long=long(:);
     386
     387    %project lat,long:
     388    [x,y]=ll2xy(lat,long,-1);
     389
     390    index=[];
     391    %  C  D
     392    %  A  B
     393    for j=1:Ny-1,
     394        for i=1:Nx-1,
     395            A=(j-1)*Nx+i;
     396            B=(j-1)*Nx+i+1;
     397            C=j*Nx+i;
     398            D=j*Nx+i+1;
     399            index(end+1,:)=[A B C];
     400            index(end+1,:)=[C B D];
     401        end
     402    end
     403
     404    %fill mesh and model:
     405    md=meshconvert(md,index,x,y);
     406    md.mesh.lat=lat;
     407    md.mesh.long=long;
     408
     409    savemodel(org,md);
     410
     411end
     412% }}}
     413% {{{ MeshGeometry:
    413414if perform(org,'MeshGeometry'),
    414        
    415         loaddata(org,'Parameters');
    416         loaddata(org,'CreateMesh');
    417         loaddata(org,'Bathymetry');
    418         loaddata(org,'IceSheetGeometry');
    419 
    420         %transfer to vertices:
    421         bathymetry=bathymetry(:);
    422         iceshelf_mask=iceshelf_mask(:);
    423         ice_mask=ice_mask(:);
    424         thickness=thickness(:);
    425         draft=draft(:);
    426 
    427         %start filling some of the fields
    428         md.geometry.bed=bathymetry;
    429         md.geometry.thickness=thickness;
    430         md.geometry.base=md.geometry.bed;
    431         pos=find(iceshelf_mask); md.geometry.base(pos)=draft(pos);
    432         md.geometry.surface=md.geometry.base+md.geometry.thickness;
    433 
    434         %nothing passes icefront:
    435         pos=find(~ice_mask);
    436         md.geometry.thickness(pos)=1;
    437         md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
    438         md.geometry.base(pos)=-di*md.geometry.thickness(pos);
    439 
    440         %level sets:
    441         md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
    442         md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
    443 
    444         pos=find(ice_mask); md.mask.ice_levelset(pos)=-1;
    445         pos=find(~iceshelf_mask & ice_mask); md.mask.groundedice_levelset(pos)=1;
    446 
    447         %identify edges of grounded ice:
    448         groundedice_levelset=md.mask.groundedice_levelset;
    449         for i=1:md.mesh.numberofelements,
    450                 m=groundedice_levelset(md.mesh.elements(i,:));
    451                 if abs(sum(m))~=3,
    452                         pos=find(m==1); md.mask.groundedice_levelset(md.mesh.elements(i,pos))=0;
    453                 end
    454         end
    455 
    456         %identify edges of ice:
    457         ice_levelset=md.mask.ice_levelset;
    458         for i=1:md.mesh.numberofelements,
    459                 m=ice_levelset(md.mesh.elements(i,:));
    460                 if abs(sum(m))~=3,
    461                         pos=find(m==-1); md.mask.ice_levelset(md.mesh.elements(i,pos))=0;
    462                 end
    463         end
    464 
    465         savemodel(org,md);
    466 end
    467 % }}}
    468 % {{{ ParameterizeIce: 
     415
     416    loaddata(org,'Parameters');
     417    loaddata(org,'CreateMesh');
     418    loaddata(org,'Bathymetry');
     419    loaddata(org,'IceSheetGeometry');
     420
     421    %transfer to vertices:
     422    bathymetry=bathymetry(:);
     423    iceshelf_mask=iceshelf_mask(:);
     424    ice_mask=ice_mask(:);
     425    thickness=thickness(:);
     426    draft=draft(:);
     427
     428    %start filling some of the fields
     429    md.geometry.bed=bathymetry;
     430    md.geometry.thickness=thickness;
     431    md.geometry.base=md.geometry.bed;
     432    pos=find(iceshelf_mask); md.geometry.base(pos)=draft(pos);
     433    md.geometry.surface=md.geometry.base+md.geometry.thickness;
     434
     435    %nothing passes icefront:
     436    pos=find(~ice_mask);
     437    md.geometry.thickness(pos)=1;
     438    md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
     439    md.geometry.base(pos)=-di*md.geometry.thickness(pos);
     440
     441    %level sets:
     442    md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
     443    md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
     444
     445    pos=find(ice_mask); md.mask.ice_levelset(pos)=-1;
     446    pos=find(~iceshelf_mask & ice_mask); md.mask.groundedice_levelset(pos)=1;
     447
     448    %identify edges of grounded ice:
     449    groundedice_levelset=md.mask.groundedice_levelset;
     450    for i=1:md.mesh.numberofelements,
     451        m=groundedice_levelset(md.mesh.elements(i,:));
     452        if abs(sum(m))~=3,
     453            pos=find(m==1); md.mask.groundedice_levelset(md.mesh.elements(i,pos))=0;
     454        end
     455    end
     456
     457    %identify edges of ice:
     458    ice_levelset=md.mask.ice_levelset;
     459    for i=1:md.mesh.numberofelements,
     460        m=ice_levelset(md.mesh.elements(i,:));
     461        if abs(sum(m))~=3,
     462            pos=find(m==-1); md.mask.ice_levelset(md.mesh.elements(i,pos))=0;
     463        end
     464    end
     465
     466    savemodel(org,md);
     467end
     468% }}}
     469% {{{ ParameterizeIce:
    469470if perform(org,'ParameterizeIce'),
    470        
    471         loaddata(org,'Parameters');
    472         loaddata(org,'CreateMesh');
    473         loaddata(org,'MeshGeometry');
    474 
    475         %miscellaneous
    476         md.miscellaneous.name='test4002';
    477 
    478         %initial velocity:
    479         md.initialization.vx=zeros(md.mesh.numberofvertices,1);
    480         md.initialization.vy=zeros(md.mesh.numberofvertices,1);
    481         md.initialization.vz=zeros(md.mesh.numberofvertices,1);
    482 
    483         %friction:
    484         md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
    485         pos=find(md.mask.groundedice_levelset<=0);
    486         md.friction.coefficient(pos)=0;
    487         md.friction.p=ones(md.mesh.numberofelements,1);
    488         md.friction.q=ones(md.mesh.numberofelements,1);
    489 
    490         %temperatures and surface mass balance:
    491         md.initialization.temperature=(273.15-20)*ones(md.mesh.numberofvertices,1);
    492         md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
    493         md.smb.mass_balance = [1*ones(md.mesh.numberofvertices,1); 1];
    494 
    495         %Flow law
    496         md.materials.rheology_B=paterson(md.initialization.temperature);
    497         md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    498         md.damage.D=zeros(md.mesh.numberofvertices,1);
    499         md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
    500        
    501         %the spcs going
    502         md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
    503         md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
    504         md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
    505         md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
    506         md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
    507         md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
    508 
    509         %deal with water:
    510         pos=find(md.mask.ice_levelset>0);
    511         md.stressbalance.spcvx(pos)=0;
    512         md.stressbalance.spcvy(pos)=0;
    513         md.stressbalance.spcvz(pos)=0;
    514         md.masstransport.spcthickness(pos)=0;
    515 
    516         %get some flux at the ice divide:
    517         pos=find(md.mesh.lat==min(md.mesh.lat));
    518         md.stressbalance.spcvy(pos)=200;
    519 
    520         %deal with boundaries, excluding icefront:
    521         vertex_on_boundary=zeros(md.mesh.numberofvertices,1);
    522         vertex_on_boundary(md.mesh.segments(:,1:2))=1;
    523         pos=find(vertex_on_boundary & md.mask.groundedice_levelset<=0);
    524         md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
    525         md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
    526         md.stressbalance.spcvz(pos)=md.initialization.vz(pos);
    527         md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
    528 
    529         md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
    530         md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
    531         md.thermal.spctemperature=[md.initialization.temperature; 1]; %impose observed temperature on surface
    532         md.basalforcings.geothermalflux=.064*ones(md.mesh.numberofvertices,1);
    533 
    534         %flow equations:
    535         md=setflowequation(md,'SSA','all');
    536 
    537         savemodel(org,md);
    538 end
    539 % }}}
    540 % {{{ RunUncoupledISSM: 
     471
     472    loaddata(org,'Parameters');
     473    loaddata(org,'CreateMesh');
     474    loaddata(org,'MeshGeometry');
     475
     476    %miscellaneous
     477    md.miscellaneous.name='test4002';
     478
     479    %initial velocity:
     480    md.initialization.vx=zeros(md.mesh.numberofvertices,1);
     481    md.initialization.vy=zeros(md.mesh.numberofvertices,1);
     482    md.initialization.vz=zeros(md.mesh.numberofvertices,1);
     483
     484    %friction:
     485    md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
     486    pos=find(md.mask.groundedice_levelset<=0);
     487    md.friction.coefficient(pos)=0;
     488    md.friction.p=ones(md.mesh.numberofelements,1);
     489    md.friction.q=ones(md.mesh.numberofelements,1);
     490
     491    %temperatures and surface mass balance:
     492    md.initialization.temperature=(273.15-20)*ones(md.mesh.numberofvertices,1);
     493    md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
     494    md.smb.mass_balance = [1*ones(md.mesh.numberofvertices,1); 1];
     495
     496    %Flow law
     497    md.materials.rheology_B=paterson(md.initialization.temperature);
     498    md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     499    md.damage.D=zeros(md.mesh.numberofvertices,1);
     500    md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
     501
     502    %the spcs going
     503    md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
     504    md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
     505    md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
     506    md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
     507    md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
     508    md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
     509
     510    %deal with water:
     511    pos=find(md.mask.ice_levelset>0);
     512    md.stressbalance.spcvx(pos)=0;
     513    md.stressbalance.spcvy(pos)=0;
     514    md.stressbalance.spcvz(pos)=0;
     515    md.masstransport.spcthickness(pos)=0;
     516
     517    %get some flux at the ice divide:
     518    pos=find(md.mesh.lat==min(md.mesh.lat));
     519    md.stressbalance.spcvy(pos)=200;
     520
     521    %deal with boundaries, excluding icefront:
     522    vertex_on_boundary=zeros(md.mesh.numberofvertices,1);
     523    vertex_on_boundary(md.mesh.segments(:,1:2))=1;
     524    pos=find(vertex_on_boundary & md.mask.groundedice_levelset<=0);
     525    md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
     526    md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
     527    md.stressbalance.spcvz(pos)=md.initialization.vz(pos);
     528    md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
     529
     530    md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     531    md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
     532    md.thermal.spctemperature=[md.initialization.temperature; 1]; %impose observed temperature on surface
     533    md.basalforcings.geothermalflux=.064*ones(md.mesh.numberofvertices,1);
     534
     535    %flow equations:
     536    md=setflowequation(md,'SSA','all');
     537
     538    savemodel(org,md);
     539end
     540% }}}
     541% {{{ RunUncoupledISSM:
    541542if perform(org,'RunUncoupledISSM'),
    542        
    543         loaddata(org,'Parameters');
    544         loaddata(org,'ParameterizeIce');
    545 
    546         %timestepping:
    547         md.timestepping.final_time=final_time;
    548         md.timestepping.time_step=time_step;
    549         md.transient.isgroundingline=1;
    550         md.transient.isthermal=0;
    551         md.groundingline.migration='SubelementMigration';
    552         md.groundingline.melt_interpolation='SubelementMelt2';
    553         md.groundingline.friction_interpolation='SubelementFriction2';
    554 
    555         md.cluster=generic('name',oshostname(),'np',2);
    556         md=solve(md,'Transient');
    557 
    558         savemodel(org,md);
     543
     544    loaddata(org,'Parameters');
     545    loaddata(org,'ParameterizeIce');
     546
     547    %timestepping:
     548    md.timestepping.final_time=final_time;
     549    md.timestepping.time_step=time_step;
     550    md.transient.isgroundingline=1;
     551    md.transient.isthermal=0;
     552    md.groundingline.migration='SubelementMigration';
     553    md.groundingline.melt_interpolation='SubelementMelt2';
     554    md.groundingline.friction_interpolation='SubelementFriction2';
     555
     556    md.cluster=generic('name',oshostname(),'np',2);
     557    md=solve(md,'Transient');
     558
     559    savemodel(org,md);
    559560end
    560561% }}}
    561562
    562563%Run MITgcm/ISSM
    563 % {{{ RunCoupledMITgcmISSM: 
     564% {{{ RunCoupledMITgcmISSM:
    564565if perform(org,'RunCoupledMITgcmISSM'),
    565566
    566         %load data:
    567         loaddata(org,'Parameters');
    568         loaddata(org,'ParameterizeIce');
    569         loaddata(org,'Bathymetry');
    570         loaddata(org,'IceSheetGeometry');
     567    %load data:
     568    loaddata(org,'Parameters');
     569    loaddata(org,'ParameterizeIce');
     570    loaddata(org,'Bathymetry');
     571    loaddata(org,'IceSheetGeometry');
    571572        endtime = round(MITgcmDeltaT * ...
    572573         floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
    573574
    574         % {{{ prepare MITgcm 
     575        % {{{ prepare MITgcm
    575576        % rename previous run directory and create new one
    576577        if exist ('run.old')
    … …  
    631632            tline = fgetl(fidi);
    632633            if ~ischar(tline), break, end
    633             %do the change here: 
     634            %do the change here:
    634635            if strcmpi(tline,' OB_Iwest = 40*1,'),
    635636                fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
    … …  
    649650        writebin('run/bathymetry.bin',bathymetry);
    650651        % }}}
    651            
    652         % {{{ ISSM settings:
    653 
    654         setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib')
    655         %timestepping:
    656         md.timestepping.start_time=start_time;
    657         md.timestepping.final_time=final_time;
    658         md.timestepping.time_step=time_step;
    659         md.cluster=generic('name',oshostname(),'np',2);
    660         md.results.TransientSolution.Base=md.geometry.base;
    661         md.transient.isgroundingline=1;
    662         md.transient.isthermal=0;
    663         md.groundingline.migration='SubelementMigration';
    664         md.groundingline.melt_interpolation='SubelementMelt2';
    665         md.groundingline.friction_interpolation='SubelementFriction2';
    666 
    667         % }}}
    668 
    669         %start looping:
    670         results=md.results;
    671 
    672         for t=start_time:time_step:final_time
     652
     653    % {{{ ISSM settings:
     654
     655    setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib')
     656    %timestepping:
     657    md.timestepping.start_time=start_time;
     658    md.timestepping.final_time=final_time;
     659    md.timestepping.time_step=time_step;
     660    md.cluster=generic('name',oshostname(),'np',2);
     661    md.results.TransientSolution.Base=md.geometry.base;
     662    md.transient.isgroundingline=1;
     663    md.transient.isthermal=0;
     664    md.groundingline.migration='SubelementMigration';
     665    md.groundingline.melt_interpolation='SubelementMelt2';
     666    md.groundingline.friction_interpolation='SubelementFriction2';
     667
     668    % }}}
     669
     670    %start looping:
     671    results=md.results;
     672
     673    for t=start_time:time_step:final_time
    673674            disp(['Year: ' num2str(t)])
    674675
    675             %send draft from ISSM to MITgcm:
    676             draft=md.results.TransientSolution(end).Base;
    677             pos=find(md.mask.ice_levelset>0); draft(pos)=0;
     676        %send draft from ISSM to MITgcm:
     677        draft=md.results.TransientSolution(end).Base;
     678        pos=find(md.mask.ice_levelset>0); draft(pos)=0;
    678679            if t>start_time
    679680                old_draft=readbin('run/icetopo.bin',[Nx,Ny]);
    … …  
    681682            writebin('run/icetopo.bin',draft);
    682683
    683             % {{{ generate MITgcm parameter file data
    684             fidi=fopen('../MITgcm/input/data','r');
    685             fido=fopen('run/data','w');
    686             tline = fgetl(fidi);
    687             fprintf(fido,'%s\n',tline);
     684        % {{{ generate MITgcm parameter file data
     685        fidi=fopen('../MITgcm/input/data','r');
     686        fido=fopen('run/data','w');
     687        tline = fgetl(fidi);
     688        fprintf(fido,'%s\n',tline);
    688689            while 1
    689690                tline = fgetl(fidi);
    690691                if ~ischar(tline), break, end
    691                 %do the change here:
    692                 if strcmpi(tline,' xgOrigin = 0.0,'),
    693                     fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
    694                     continue;
    695                 end
    696                 if strcmpi(tline,' ygOrigin = -80.0,'),
    697                     fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
    698                     continue;
    699                 end
    700                 if strcmpi(tline,' delX = 20*0.25,'),
    701                     fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
    702                     continue;
    703                 end
    704                 if strcmpi(tline,' delY = 20*0.25,'),
    705                     fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
    706                     continue;
    707                 end
     692        %do the change here:
     693        if strcmpi(tline,' xgOrigin = 0.0,'),
     694            fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
     695            continue;
     696        end
     697        if strcmpi(tline,' ygOrigin = -80.0,'),
     698            fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
     699            continue;
     700        end
     701        if strcmpi(tline,' delX = 20*0.25,'),
     702            fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
     703            continue;
     704        end
     705        if strcmpi(tline,' delY = 20*0.25,'),
     706            fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
     707            continue;
     708        end
    708709                if strcmpi(tline,' delZ = 30*30.0,'),
    709710                    fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
    710711                    continue;
    711712                end
    712                 if strcmpi(tline,' endTime=2592000.,'),
    713                     fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
    714                     continue;
    715                 end
    716                 if strcmpi(tline,' deltaT=1200.0,'),
    717                     fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
    718                     continue;
    719                 end
    720                 if strcmpi(tline,' pChkptFreq=2592000.,'),
    721                     fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
    722                     continue;
    723                 end
    724                 if strcmpi(tline,' taveFreq=2592000.,'),
    725                     fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
    726                     continue;
    727                 end
     713        if strcmpi(tline,' endTime=2592000.,'),
     714            fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
     715            continue;
     716        end
     717        if strcmpi(tline,' deltaT=1200.0,'),
     718            fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
     719            continue;
     720        end
     721        if strcmpi(tline,' pChkptFreq=2592000.,'),
     722            fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
     723            continue;
     724        end
     725        if strcmpi(tline,' taveFreq=2592000.,'),
     726            fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
     727            continue;
     728        end
    728729                if strcmpi(tline,' rhoConst=1030.,'),
    729730                    fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
    … …  
    734735                    continue;
    735736                end
    736                 fprintf(fido,'%s\n',tline);
    737             end
    738             %close  files
    739             fclose(fidi);
    740             fclose(fido);
    741             % }}}
    742 
    743             % {{{ generate initial MITgcm conditions
     737        fprintf(fido,'%s\n',tline);
     738        end
     739        %close  files
     740        fclose(fidi);
     741        fclose(fido);
     742        % }}}
     743
     744        % {{{ generate initial MITgcm conditions
    744745            ds=round(endtime/MITgcmDeltaT);
    745746            if t>start_time
    … …  
    760761                tmp(find(tmp<0))=0;
    761762                [im jm]=find(tmp); % horizontal indices where there is melt
    762                
     763
    763764                % Extrapolate T/S to locations where ice shelf retreated
    764765                for i=1:length(im)
    … …  
    794795            % }}}
    795796
    796             % {{{ system call to run MITgcm 
     797            % {{{ system call to run MITgcm
    797798            cd run
    798799            eval(['!mpiexec -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
    … …  
    810811            % }}}
    811812
    812             %get melting rates from MITgcm
    813             %upward fresh water flux (kg/m^2/s):
    814             fnm=['run/SHICE_fwFluxtave_' myint2str(ts,10) '.data'];
    815             melting_rate=readbin(fnm,[Nx Ny]);
    816 
    817             %send averaged melting rate to ISSM
    818             %downward fresh water flux (m/y):
    819             melting_rate=-melting_rate(:)*y2s/rho_ice;
    820                 md.basalforcings.floatingice_melting_rate=melting_rate;
    821 
    822             % {{{ run ISSM and recover results
    823 
    824             md.timestepping.start_time=t;
    825             md.timestepping.final_time=t+time_step;;
    826                 md=solve(md,'Transient');
    827 
    828                 base=md.results.TransientSolution(end).Base;
    829                 thickness=md.results.TransientSolution(end).Thickness;
    830                 md.geometry.base=base;
    831                 md.geometry.thickness=thickness;
    832                 md.geometry.surface=md.geometry.base+md.geometry.thickness;
    833                 md.initialization.vx=md.results.TransientSolution(end).Vx;
    834                 md.initialization.vy=md.results.TransientSolution(end).Vy;
    835                 md.initialization.vel=md.results.TransientSolution(end).Vel;
    836                 md.initialization.pressure=md.results.TransientSolution(end).Pressure;
    837                 md.mask.groundedice_levelset=md.results.TransientSolution(end).MaskGroundediceLevelset;
    838                 md.results.TransientSolution(end).FloatingiceMeltingRate=md.basalforcings.floatingice_melting_rate;
    839                
    840                 %save these results in the model, otherwise, they'll be wiped out
    841                 results(end+1)=md.results;
    842 
    843                 % }}}
    844 
    845 
    846         end
    847 
    848         md.results=results;
    849         savemodel(org,md);
    850 end
    851 % }}}
    852 % {{{ RunCoupledMITgcmISSM2: 
     813        %get melting rates from MITgcm
     814        %upward fresh water flux (kg/m^2/s):
     815        fnm=['run/SHICE_fwFluxtave_' myint2str(ts,10) '.data'];
     816        melting_rate=readbin(fnm,[Nx Ny]);
     817
     818        %send averaged melting rate to ISSM
     819        %downward fresh water flux (m/y):
     820        melting_rate=-melting_rate(:)*y2s/rho_ice;
     821        md.basalforcings.floatingice_melting_rate=melting_rate;
     822
     823        % {{{ run ISSM and recover results
     824
     825        md.timestepping.start_time=t;
     826        md.timestepping.final_time=t+time_step;;
     827        md=solve(md,'Transient');
     828
     829        base=md.results.TransientSolution(end).Base;
     830        thickness=md.results.TransientSolution(end).Thickness;
     831        md.geometry.base=base;
     832        md.geometry.thickness=thickness;
     833        md.geometry.surface=md.geometry.base+md.geometry.thickness;
     834        md.initialization.vx=md.results.TransientSolution(end).Vx;
     835        md.initialization.vy=md.results.TransientSolution(end).Vy;
     836        md.initialization.vel=md.results.TransientSolution(end).Vel;
     837        md.initialization.pressure=md.results.TransientSolution(end).Pressure;
     838        md.mask.groundedice_levelset=md.results.TransientSolution(end).MaskGroundediceLevelset;
     839        md.results.TransientSolution(end).FloatingiceMeltingRate=md.basalforcings.floatingice_melting_rate;
     840
     841        %save these results in the model, otherwise, they'll be wiped out
     842        results(end+1)=md.results;
     843
     844        % }}}
     845
     846
     847    end
     848
     849    md.results=results;
     850    savemodel(org,md);
     851end
     852% }}}
     853% {{{ RunCoupledMITgcmISSM2:
    853854if perform(org,'RunCoupledMITgcmISSM2'),
    854        
    855         loaddata(org,'Parameters');
    856         loaddata(org,'ParameterizeIce');
    857         loaddata(org,'Bathymetry');
    858         loaddata(org,'IceSheetGeometry');
    859                 endtime = round(MITgcmDeltaT * floor(final_time*y2s/MITgcmDeltaT));
    860                 outputtime = round(MITgcmDeltaT * floor(time_step*y2s/MITgcmDeltaT));
    861 
    862         % {{{ prepare MITgcm 
     855
     856    loaddata(org,'Parameters');
     857    loaddata(org,'ParameterizeIce');
     858    loaddata(org,'Bathymetry');
     859    loaddata(org,'IceSheetGeometry');
     860        endtime = round(MITgcmDeltaT * floor(final_time*y2s/MITgcmDeltaT));
     861        outputtime = round(MITgcmDeltaT * floor(time_step*y2s/MITgcmDeltaT));
     862
     863        % {{{ prepare MITgcm
    863864        % rename previous run directory and create new one
    864865        if exist ('run.old')
    … …  
    918919            tline = fgetl(fidi);
    919920            if ~ischar(tline), break, end
    920             %do the change here: 
     921            %do the change here:
    921922            if strcmpi(tline,' OB_Iwest = 40*1,'),
    922923                fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
    … …  
    937938        writebin('run/icetopo.bin',draft);
    938939        % }}}
    939         % {{{ generate MITgcm parameter file data 
     940        % {{{ generate MITgcm parameter file data
    940941        fidi=fopen('../MITgcm/input/data','r');
    941942        fido=fopen('run/data','w');
    … …  
    945946            tline = fgetl(fidi);
    946947            if ~ischar(tline), break, end
    947             %do the change here: 
     948            %do the change here:
    948949            if strcmpi(tline,' xgOrigin = 0.0,'),
    949950                fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
    … …  
    996997        fclose(fido);
    997998        % }}}
    998                
    999         md.transient.isoceancoupling=1;
    1000         md.transient.isgroundingline=1;
    1001         md.groundingline.migration='None';
    1002         md.groundingline.melt_interpolation='SubelementMelt2';
    1003         md.groundingline.friction_interpolation='SubelementFriction2';
    1004         md.timestepping.coupling_time=time_step;
     999
     1000    md.transient.isoceancoupling=1;
     1001    md.transient.isgroundingline=1;
     1002    md.groundingline.migration='None';
     1003    md.groundingline.melt_interpolation='SubelementMelt2';
     1004    md.groundingline.friction_interpolation='SubelementFriction2';
     1005    md.timestepping.coupling_time=time_step;
    10051006   md.timestepping.time_step=time_step;
    1006         md.timestepping.final_time=final_time-time_step;
    1007         md.cluster.npocean=nPx*nPy;
    1008         md.cluster.np=2;
    1009         md.cluster.executionpath=[pwd '/run'];
    1010         md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};
    1011 
    1012         md=solveiceocean(md,'Transient','runtimename',false);
    1013 
    1014 %       %eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' ']);
    1015 %       eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' : -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
     1007    md.timestepping.final_time=final_time-time_step;
     1008    md.cluster.npocean=nPx*nPy;
     1009    md.cluster.np=2;
     1010    md.cluster.executionpath=[pwd '/run'];
     1011    md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};
     1012
     1013    md=solveiceocean(md,'Transient','runtimename',false);
     1014
     1015%   %eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' ']);
     1016%   eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' : -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
    10161017end
    10171018% }}}
    … …  
    10271028melting_rate_4=readbin(fnm,[Nx Ny]);
    10281029field_names     ={'Base1','Melting1','Vx2','Vy2','Thickness2','Base2','MaskGroundediceLevelset2','FloatingiceMeltingRate2',...
    1029         'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',...
    1030         'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'};
     1030    'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',...
     1031    'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'};
    10311032field_tolerances={2e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,...
    1032         1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,...
    1033         1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 };
     1033    1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13,...
     1034    1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 };
    10341035field_values={...
    1035         (md.results.TransientSolution(1).Base),...
    1036         (melting_rate_1(:)),...
    1037         (md.results.TransientSolution(2).Vx),...
    1038         (md.results.TransientSolution(2).Vy),...
    1039         (md.results.TransientSolution(2).Thickness),...
    1040         (md.results.TransientSolution(2).Base),...
    1041         (md.results.TransientSolution(2).MaskGroundediceLevelset),...
    1042         (md.results.TransientSolution(2).BasalforcingsFloatingiceMeltingRate),...
    1043         (melting_rate_2(:)),...
    1044         (md.results.TransientSolution(3).Vx),...
    1045         (md.results.TransientSolution(3).Vy),...
    1046         (md.results.TransientSolution(3).Thickness),...
    1047         (md.results.TransientSolution(3).Base),...
    1048         (md.results.TransientSolution(3).MaskGroundediceLevelset),...
    1049         (md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),...
    1050         (melting_rate_3(:)),...
    1051         (md.results.TransientSolution(4).Vx),...
    1052         (md.results.TransientSolution(4).Vy),...
    1053         (md.results.TransientSolution(4).Thickness),...
    1054         (md.results.TransientSolution(4).Base),...
    1055         (md.results.TransientSolution(4).MaskGroundediceLevelset),...
    1056         (md.results.TransientSolution(4).BasalforcingsFloatingiceMeltingRate),...
    1057         (melting_rate_4(:)),...
    1058         };
     1036    (md.results.TransientSolution(1).Base),...
     1037    (melting_rate_1(:)),...
     1038    (md.results.TransientSolution(2).Vx),...
     1039    (md.results.TransientSolution(2).Vy),...
     1040    (md.results.TransientSolution(2).Thickness),...
     1041    (md.results.TransientSolution(2).Base),...
     1042    (md.results.TransientSolution(2).MaskGroundediceLevelset),...
     1043    (md.results.TransientSolution(2).BasalforcingsFloatingiceMeltingRate),...
     1044    (melting_rate_2(:)),...
     1045    (md.results.TransientSolution(3).Vx),...
     1046    (md.results.TransientSolution(3).Vy),...
     1047    (md.results.TransientSolution(3).Thickness),...
     1048    (md.results.TransientSolution(3).Base),...
     1049    (md.results.TransientSolution(3).MaskGroundediceLevelset),...
     1050    (md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),...
     1051    (melting_rate_3(:)),...
     1052    (md.results.TransientSolution(4).Vx),...
     1053    (md.results.TransientSolution(4).Vy),...
     1054    (md.results.TransientSolution(4).Thickness),...
     1055    (md.results.TransientSolution(4).Base),...
     1056    (md.results.TransientSolution(4).MaskGroundediceLevelset),...
     1057    (md.results.TransientSolution(4).BasalforcingsFloatingiceMeltingRate),...
     1058    (melting_rate_4(:)),...
     1059    };

  • issm/trunk-jpl/test/NightlyRun/test4003.m

    r24257 r24456  
    33%
    44%Script control parameters
    5 steps=[1 2 3 4 5 6 7 8 9 10 11];
     5steps=1:11;
    66final_time=1/365;
    77
    88%To download and recompile MITgcm from scratch:
    9 !rm -rf $ISSM_DIR/test/MITgcm/install_dngoldberg
    10 !rm -rf $ISSM_DIR/test/MITgcm/build/*
     9!rm -rf ${ISSM_DIR}/test/MITgcm/install_dngoldberg
     10!rm -rf ${ISSM_DIR}/test/MITgcm/build/*
     11!rm -rf Models
    1112
    1213%Organizer
    13 mkdir Models
    14 org=organizer('repository','Models/','prefix','IceOcean.','steps',steps);
     14!mkdir Models
     15org=organizer('repository','Models','prefix','IceOcean.','steps',steps);
    1516
    1617presentdirectory=pwd;
Note: See TracChangeset for help on using the changeset viewer.