Changeset 26719


Ignore:
Timestamp:
12/08/21 17:08:48 (3 years ago)
Author:
dmenemen
Message:

updated test4003 with 3*200 domain based on
MITgcm/verification/shelfice_2d_remesh
runs but hangs on step 11: RunCoupledMITgcmISSM

Location:
issm/trunk-jpl/test
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/test/MITgcm/code_4003/cpl_issm.F

    r24121 r26719  
    6767     &        MPI_COMM_MODEL,mpiRC)
    6868         _END_MASTER( myThid )
    69 C        print*, 'Ocean received CouplingTime: ', CouplingTime
     69#ifdef ALLOW_DEBUG
     70         IF ( debugLevel .GE. debLevD ) THEN
     71            print*, 'Ocean received CouplingTime: ', CouplingTime
     72         ENDIF
     73#endif
    7074
    7175C   I2. MITgcm sends grid size (NX and NY)
    … …  
    135139            call MPI_Recv(IceModelTime,1,MPI_DOUBLE,0,10001001,
    136140     &           toissmcomm,mpistatus,mpiRC)
    137 C           print*, 'Ocean received IceModelTime: ', IceModelTime
     141#ifdef ALLOW_DEBUG
     142            IF ( debugLevel .GE. debLevD ) THEN
     143               print*, 'Ocean received IceModelTime: ', IceModelTime
     144            ENDIF
     145#endif
    138146            _END_MASTER( myThid )
    139147         ENDIF
    … …  
    176184         ENDIF
    177185         CALL BAR2( myThid )
    178 C        print*,'Done Sending shelficeFreshWaterFlux array.'
    179          
     186#ifdef ALLOW_DEBUG
     187         IF ( debugLevel .GE. debLevD ) THEN
     188            print*,'Done Sending shelficeFreshWaterFlux array.'
     189         ENDIF
     190#endif
    180191      ENDIF
    181192C End recurring step C3.

  • issm/trunk-jpl/test/MITgcm/input_4003/data

    r26664 r26719  
    4646#-
    4747 readBinaryPrec=64,
    48  writeBinaryPrec=64,
     48#writeBinaryPrec=64,
    4949 useSingleCpuIO=.TRUE.,
    5050#globalFiles=.TRUE.,

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

    r26664 r26719  
    33%
    44%Script control parameters
    5 steps=1:12;
    6 steps=1:6;
    7 final_time=1/365;
     5steps=1:11;
    86
    97%To download and recompile MITgcm from scratch:
    … …  
    2018% {{{ Parameters:
    2119if perform(org,'Parameters'),
    22         Nx=3; %number of longitude cells
    23         Ny=200; %number of latitude cells
    24         Nz=90; %number of MITgcm vertical cells
    25         nPx=1; %number of MITgcm processes to use in x direction
    26         nPy=8; %number of MITgcm processes to use in y direction
    27         xgOrigin=0; %origin of longitude
    28         ygOrigin=-75.5; %origin of latitude
    29         dLong=.125; %longitude grid spacing
    30         dLat=dLong/16; %latitude grid spacing
    31         delZ=10; %thickness of vertical levels
    32         icefront_position_ratio=.75;
    33         ice_thickness=1000;
    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=450; % 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         % 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);
     20    Nx=3;   % number of longitude cells
     21    Ny=200; % number of latitude cells
     22    Nz=90;  % number of MITgcm vertical cells
     23    nPx=1;  % number of MITgcm processes to use in x direction
     24    nPy=8;  % number of MITgcm processes to use in y direction
     25    xgOrigin=0;     % origin of longitude
     26    ygOrigin=-75.5; % origin of latitude
     27    dLong=.125;     % longitude grid spacing
     28    dLat=dLong/16;  % latitude grid spacing
     29    delZ=10;        % thickness of vertical levels (m)
     30    gravity= 9.81;  % gravity (m^2/s)
     31    rho_ice=917;
     32    rho_water=1030;
     33    di=rho_ice/rho_water;
     34    prec = 'real*8'; % precision of MITgcm input binary files
     35
     36    % bathymetry and ice sheet geometry
     37    H = -900;       % water depth in the ice shelf cavity
     38    Hmin = -600;    % deepest point of cavern           
     39    Hmax = -300;    % shallowest point of cavern
     40    jEnd = Ny*3/4;  % where ice-shelf ends
     41    j2   = jEnd+1;
     42
     43    savedata(org, Nx, Ny, nPx, nPy, Nz, dLong, dLat, delZ, xgOrigin, ygOrigin, ...
     44             rho_ice, rho_water, di, H, Hmin, Hmax, jEnd, j2, gravity, prec);
    7145end
    7246% }}}
    … …  
    7549
    7650    loaddata(org,'Parameters');
    77 
    78     v0 = 2e3;
    79     h0 = 800;
    80 
    81     hfacMin = 0.2;
    82 
    83     eos = 'jmd95z';
    84     prec = 'real*8';
    8551
    8652    %create lat,lon
    … …  
    8955    long = xgOrigin+[0:Nx-1]*dLong;
    9056    lonc = long+dLong/2;
     57    [lat lon]=meshgrid(latc,lonc);
    9158    zC=-delZ*([1:Nz]-0.5);
    9259    zF=-delZ*[0:Nz];
    93            
    94     lat=(ygOrigin+dLat/2):dLat:(ygOrigin+Ny*dLat);
    95     lon=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong);
    96     [lat lon]=meshgrid(lat,lon);
    97 
    98     longmin=min(lon(:));
    99     longmax=max(lon(:));
    100     latmin=min(lat(:));
    101     latmax=max(lat(:));
    102 
    103 % Gravity
    104 gravity= 9.81;
    105 rhoConst= 1030;
    106 
    107 % Nominal depth of model (meters)
    108 H = -900;               %water depth in the ice shelf cavity
    109 Hmin = -600;            % deepest point of cavern               
    110 Hmax = -300;            % shallowest point of cavern
    111 jEnd = Ny*3/4;           % where ice-shelf ends
    112 dHdy = (Hmax-Hmin)/dLat/(jEnd-2); %Slope of ice shelf
    113 
     60   
    11461    %create bathymetry:
    11562    bathymetry = ones(Nx,Ny)*H;
    … …  
    11966    savedata(org,lat,lon,bathymetry);
    12067
    121 j2=jEnd+1;
    122 hIce=bathymetry;
    123 for i=1:Nx
    124     hIce(i,:)=Hmin+dHdy*[-1:Ny-2]*dLat;
    125     hIce(i,1)=0; hIce(i,j2:Ny)=0;
    126 end
    127 
    128 var=([1:Ny]-2)/(jEnd-2);
    129 dMdt_fy=bathymetry;
    130 for i=1:Nx
    131     dMdt_fy(i,:)=-cos(pi*var);
    132     dMdt_fy(i,1)=0; dMdt_fy(i,j2:Ny)=0;
    133 end
    134 
    135 regMsk=ones(Nx,Ny);
    136 regMsk(:,1)=0; regMsk(:,j2:Ny)=2;
    137 
    138 %- rate of change due to ice-stream dynamics
    139 rateDyn=rhoConst*0.1/3600;  sfx='r02';
    140 
    141 dMdt=rateDyn*dMdt_fy;
    142 
    143 dz = delZ*ones(1,Nz);
    144 zgp1 = [0,cumsum(dz)];
    145 zc = .5*(zgp1(1:end-1)+zgp1(2:end));
    146 zg = zgp1(1:end-1);
    147 dz = diff(zgp1);
    148 
    149 T_sfc = -1.9;
    150 T_bot = 2;
    151 del_T = (T_bot - T_sfc)/(59*delZ);
    152 tref=0*dz;
    153 for k = 1:Nz;
    154     tref(k) = T_sfc + del_T*((k-20)*delZ);
    155     tref(k)= max(T_sfc,min(tref(k),T_bot));
    156 end
    157 
    158 S_sfc = 34.2;
    159 S_bot = 34.7;
    160 del_S = (S_bot - S_sfc)/(59*delZ);
    161 sref=0*dz;
    162 for k = 1:Nz;
    163     sref(k) = S_sfc + del_S*((k-20)*delZ);
    164     sref(k)= max(S_sfc,min(sref(k),S_bot));
    165 end
    166 pEOS=-rhoConst*gravity*zC; % in Pa
    167 pEOS=pEOS*1.e-4; % in dBar
    168 rhoAn=densjmd95(sref,tref,pEOS);
    169 rhoAn=rhoAn-rhoConst;
    170 
    171 pF=-rhoConst*gravity*zF*1.e-4; % in dBar
    172 rhoUp=densjmd95(sref,tref,pF(2:end));
    173 rhoDw=densjmd95(sref,tref,pF(1:Nz));
    174 dRho=rhoUp(1:Nz-1)-rhoDw(2:Nz);
    175 NSq=-gravity*dRho/delZ/rhoConst;
    176 
    177 mnV=min(NSq); MxV=max(NSq); Avr=mean(NSq);
    178 
    179 zax=[1:Nz];
    180 
    181 v1=2.5e-2;
    182 var=1+Nz-2*zax; var=var/(Nz-1);
    183 vobc=v1*var;
    184 
    185 rhoAvr=rhoConst-1.345;
    186 mIce=-rhoAvr*hIce;
    187 
    18868end
    18969% }}}
    … …  
    19171if perform(org,'IceSheetGeometry'),
    19272
    193         loaddata(org,'Parameters');
    194         loaddata(org,'Bathymetry');
    195         latmin=min(lat(:));
    196         latmax=max(lat(:));
    197 
    198         %put ice_thickness constant layer of ice over the bathymetry, unless it floats:
    199         s=size(bathymetry);
    200         thickness=ice_thickness*ones(s);
    201 
    202         %figure out ice shelf:
    203         pos=find(-di*thickness>bathymetry);
    204         iceshelf_mask=zeros(s);
    205         iceshelf_mask(pos)=1;
    206 
    207         ice_mask=ones(s);
    208         pos=find((lat-latmin)/(latmax-latmin)>(icefront_position_ratio));
    209         ice_mask(pos)=0;
    210         iceshelf_mask(pos)=0;
    211 
    212         %compute draft of ice shelf:
    213         draft=bathymetry;
    214         pos=find(iceshelf_mask);
    215         draft(pos)=-di*thickness(pos);
    216         pos=find(~ice_mask);
    217         draft(pos)=0;
    218 
    219         savedata(org,ice_mask,iceshelf_mask,draft,thickness);
     73    loaddata(org,'Parameters');
     74    loaddata(org,'Bathymetry');
     75    latmin=min(lat(:));
     76    latmax=max(lat(:));
     77
     78    dHdy = (Hmax-Hmin)/dLat/(jEnd-2); %Slope of ice shelf
     79    draft=bathymetry;
     80    for i=1:Nx
     81        draft(i,:)=Hmin+dHdy*[-1:Ny-2]*dLat;
     82        draft(i,1)=0; draft(i,j2:Ny)=0;
     83    end
     84
     85    ice_mask=ones(Nx,Ny);
     86    ice_mask(:,j2:Ny)=0;
     87    iceshelf_mask=ice_mask;
     88    iceshelf_mask(:,1)=0;
     89    thickness=abs(draft)/di;
     90    thickness(:,1)=thickness(:,2)+draft(:,2);
     91   
     92    savedata(org,ice_mask,iceshelf_mask,draft,thickness);
    22093end
    22194% }}}
    … …  
    235108% {{{ RunUncoupledMITgcm:
    236109if perform(org,'RunUncoupledMITgcm'),
     110
     111    loaddata(org,'Parameters');
     112    loaddata(org,'Bathymetry');
     113    loaddata(org,'IceSheetGeometry');
     114
    237115    % rename previous run directory and create new one
    238116    if exist ('run.old')
    … …  
    246124    !\cp ../MITgcm/input_4003/* run
    247125    !\cp ../MITgcm/input_4003/eedata_uncoupled run/eedata
    248     cd run
    249    
    250     namF='bathy_flat.bin';
     126   
     127    % {{{ Construct MITgcm binary input files
     128    namF='run/bathy_flat.bin';
    251129    fid=fopen(namF,'w','b'); fwrite(fid,bathymetry,prec);fclose(fid);
    252130   
    253     namF='shelficeTopo.Lin.bin';
    254     fid=fopen(namF,'w','b'); fwrite(fid,hIce,prec);fclose(fid);
    255    
    256     namF='under_Ice_mask.bin';
    257     fid=fopen(namF,'w','b'); fwrite(fid,regMsk,prec);fclose(fid);
    258    
    259     namF=sprintf('%s.%s.%s','shelfice_dMdt',sfx,'bin');
     131    namF='run/shelficeTopo.Lin.bin';
     132    fid=fopen(namF,'w','b'); fwrite(fid,draft,prec);fclose(fid);
     133   
     134    var=([1:Ny]-2)/(jEnd-2);
     135    dMdt_fy=zeros(Nx,Ny);
     136    for i=1:Nx
     137        dMdt_fy(i,:)=-cos(pi*var);
     138        dMdt_fy(i,1)=0; dMdt_fy(i,j2:Ny)=0;
     139    end   
     140    %- rate of change due to ice-stream dynamics
     141    rateDyn=rho_water*0.1/3600;  sfx='r02';
     142    dMdt=rateDyn*dMdt_fy;
     143    namF=sprintf('%s.%s.%s','run/shelfice_dMdt',sfx,'bin');
    260144    fid=fopen(namF,'w','b'); fwrite(fid,dMdt,prec);fclose(fid);
    261145
    262     namF='temp_obc.bin';
     146    T_sfc = -1.9;
     147    T_bot = 2;
     148    del_T = (T_bot - T_sfc)/(59*delZ);
     149    tref=zeros(1,Nz);
     150    for k = 1:Nz;
     151        tref(k) = T_sfc + del_T*((k-20)*delZ);
     152        tref(k)= max(T_sfc,min(tref(k),T_bot));
     153    end
     154    namF='run/temp_obc.bin';
    263155    tref=[tref; tref; tref];
    264156    fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid);
    265157
    266     namF='salt_obc.bin';
     158    S_sfc = 34.2;
     159    S_bot = 34.7;
     160    del_S = (S_bot - S_sfc)/(59*delZ);
     161    sref=zeros(1,Nz);
     162    for k = 1:Nz;
     163        sref(k) = S_sfc + del_S*((k-20)*delZ);
     164        sref(k)= max(S_sfc,min(sref(k),S_bot));
     165    end
     166    namF='run/salt_obc.bin';
    267167    sref=[sref; sref; sref];
    268168    fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid);
    269169
    270     namF='vVel_obc.bin';
     170    zax=[1:Nz];
     171    v1=2.5e-2;
     172    var=1+Nz-2*zax; var=var/(Nz-1);
     173    vobc=v1*var;
     174    namF='run/vVel_obc.bin';
    271175    vobc=[vobc; vobc; vobc];
    272176    fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid);
    … …  
    276180            var(i,j,:)=tref(1,:);
    277181        end, end
    278     namF='temp_ini.bin';
     182    namF='run/temp_ini.bin';
    279183    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
    280184
    … …  
    282186            var(i,j,:)=sref(1,:);
    283187        end, end
    284     namF='salt_ini.bin';
     188    namF='run/salt_ini.bin';
    285189    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
    286190
    287     namF='shelficeMass.Lin.bin';
     191    rhoAvr=rho_water-1.345;
     192    mIce=-rhoAvr*draft;
     193    namF='run/shelficeMass.Lin.bin';
    288194    fid=fopen(namF,'w','b'); fwrite(fid,mIce,prec);fclose(fid);
    289 
     195    % }}}
     196
     197    cd run
    290198    eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
    291199    cd ..
    … …  
    297205if perform(org,'CreateMesh'),
    298206
    299         %create model:
    300         md=model();
    301 
    302         %Grab lat,long from MITgcm:
    303         long=readbin('run/XG.data',[Nx Ny]);
    304         long=[long long(:,end)]; long=[long; -105.1250*ones(1,size(long,2))];
    305         lat=readbin('run/YG.data',[Nx Ny]);
    306         lat=[lat -73.8832*ones(size(lat,1),1)]; lat=[lat; lat(end,:)];
    307 
    308         %project lat,long:
    309         [x,y]=ll2xy(lat(:),long(:),-1);
    310 
    311         Nx=size(lat,1); Ny=size(lat,2);
    312         index=[];
    313         %  C  D
    314         %  A  B
    315         for j=1:Ny-1,
    316                 for i=1:Nx-1,
    317                         A=(j-1)*Nx+i;
    318                         B=(j-1)*Nx+i+1;
    319                         C=j*Nx+i;
    320                         D=j*Nx+i+1;
    321                         index(end+1,:)=[A B C];
    322                         index(end+1,:)=[C B D];
    323                 end
    324         end
    325 
    326         %fill mesh and model:
    327         md=meshconvert(md,index,x,y);
    328         md.mesh.lat=lat(:);
    329         md.mesh.long=long(:);
    330 
    331         savemodel(org,md);
     207    loaddata(org,'Parameters');
     208
     209    %create model:
     210    md=model();
     211
     212    %Grab lat,long from MITgcm:
     213    long=readbin('run/XG.data',[Nx Ny]);
     214    long=[long long(:,end)]; long=[long; 2*long(Nx,:)-long(Nx-1,:)];
     215    lat=readbin('run/YG.data',[Nx Ny]);
     216    lat=[lat 2*lat(:,Ny)-lat(:,Ny-1)]; lat=[lat; lat(end,:)];
     217
     218    %project lat,long:
     219    [x,y]=ll2xy(lat(:),long(:),-1);
     220    index=[];
     221    %  C  D
     222    %  A  B
     223    for j=1:Ny,
     224        for i=1:Nx,
     225            A=(j-1)*(Nx+1)+i;
     226            B=(j-1)*(Nx+1)+i+1;
     227            C=j*(Nx+1)+i;
     228            D=j*(Nx+1)+i+1;
     229            index(end+1,:)=[A B C];
     230            index(end+1,:)=[C B D];
     231        end
     232    end
     233
     234    %fill mesh and model:
     235    md=meshconvert(md,index,x,y);
     236    md.mesh.lat=lat(:);
     237    md.mesh.long=long(:);
     238
     239    savemodel(org,md);
    332240
    333241end
    … …  
    336244if perform(org,'MeshGeometry'),
    337245
    338         loaddata(org,'CreateMesh');
    339 
    340         %transfer to vertices:
    341         bathymetry=readbin('run/bathy.box',[3 200],1,'real*8');
    342         bathymetry=[bathymetry bathymetry(:,end)]; bathymetry=[bathymetry(1,:); bathymetry];
    343         iceshelf_mask=-1*ones(size(bathymetry));
    344         ice_mask=readbin('run/hmask3.box',[3 200],1,'real*8');
    345         ice_mask=[ice_mask ice_mask(:,end)]; ice_mask=[ice_mask(1,:); ice_mask];
    346         thickness=readbin('run/h0.bin',[3 200],1,'real*8');
    347         thickness=[thickness thickness(:,end)]; thickness=[thickness; thickness(end,:)];
    348 
    349         %start filling some of the fields
    350         md.geometry.bed=bathymetry(:);
    351         md.geometry.thickness=thickness(:);
    352         md.geometry.base=-917/1028*md.geometry.thickness;
    353         md.geometry.surface=md.geometry.base+md.geometry.thickness;
    354 
    355         %nothing passes icefront:
    356         pos=find((~ice_mask(:) & ice_mask(:)~=0) | thickness(:)==0);
    357         md.geometry.thickness(pos)=1;
    358         md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
    359         md.geometry.base(pos)=-di*md.geometry.thickness(pos);
    360 
    361         %level sets:
    362         md.mask.ocean_levelset=iceshelf_mask(:);
    363         md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
    364 
    365         pos=find(~ice_mask(:) & thickness(:)==0); md.mask.ice_levelset(pos)=1;
    366 
    367         savemodel(org,md);
     246    loaddata(org,'Parameters');
     247    loaddata(org,'CreateMesh');
     248    loaddata(org,'Bathymetry');
     249    loaddata(org,'IceSheetGeometry');
     250
     251    %transfer to vertices:
     252    bathymetry=[bathymetry bathymetry(:,end)]; bathymetry=[bathymetry(1,:); bathymetry];
     253    ice_mask=[ice_mask ice_mask(:,end)]; ice_mask=[ice_mask(1,:); ice_mask];
     254    iceshelf_mask=[iceshelf_mask iceshelf_mask(:,end)]; iceshelf_mask=[iceshelf_mask(1,:); iceshelf_mask];
     255    thickness=[thickness thickness(:,end)]; thickness=[thickness; thickness(end,:)];
     256
     257    %start filling some of the fields
     258    md.geometry.bed=bathymetry(:);
     259    md.geometry.thickness=thickness(:);
     260    md.geometry.base=-917/1028*md.geometry.thickness;
     261    md.geometry.surface=md.geometry.base+md.geometry.thickness;
     262
     263    %nothing passes icefront:
     264    pos=find((~ice_mask(:) & ice_mask(:)~=0) | thickness(:)==0);
     265    md.geometry.thickness(pos)=1;
     266    md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
     267    md.geometry.base(pos)=-di*md.geometry.thickness(pos);
     268
     269    %level sets:
     270    md.mask.ocean_levelset=iceshelf_mask(:);
     271    md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
     272
     273    pos=find(~ice_mask(:) & thickness(:)==0); md.mask.ice_levelset(pos)=1;
     274
     275    savemodel(org,md);
     276
    368277end
    369278% }}}
    … …  
    463372if perform(org,'RunCoupledMITgcmISSM'),
    464373
    465         loaddata(org,'Parameters');
    466         loaddata(org,'ParameterizeIce');
    467 
    468         md=loadmodel(org,'RunUncoupledISSM');
    469 
    470                 endtime = round(MITgcmDeltaT * floor(final_time*y2s/MITgcmDeltaT));
    471                 outputtime = round(MITgcmDeltaT * floor(time_step*y2s/MITgcmDeltaT));
    472 
    473         % {{{ prepare ISSM: start from the steady-state
    474 
    475                   md.geometry.base=md.results.TransientSolution(end).Base;
    476                   md.geometry.surface=md.results.TransientSolution(end).Surface;
    477                   md.geometry.thickness=md.results.TransientSolution(end).Thickness;
    478                   md.initialization.vx=md.results.TransientSolution(end).Vx;
    479                   md.initialization.vy=md.results.TransientSolution(end).Vy;
    480                   md.initialization.vel=md.results.TransientSolution(end).Vel;
    481                   md.initialization.pressure=md.results.TransientSolution(end).Pressure;
    482                   md.transient.isoceancoupling=1;
    483                   md.transient.isgroundingline=0;
    484                   md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};
    485 
    486         % }}}
    487         % {{{ prepare MITgcm
    488         % rename previous run directory and create new one
    489         if exist ('run.old')
    490             !\rm -rf run.old
    491         end
    492         if exist ('run')
    493             !\mv run run.old
    494         end
    495         !\mkdir run
    496         !\cp ../MITgcm/build/mitgcmuv run
    497         !\cp ../MITgcm/install/verification/shelfice_remeshing/input/* run
    498         !\cp ../MITgcm/input_4003/* run
    499         % }}}
    500                   md.timestepping.coupling_time=time_step;
    501                   md.timestepping.time_step=time_step;
    502                   md.timestepping.final_time=final_time-time_step;
    503                   md.cluster.npocean=nPx*nPy;
    504                   md.cluster.np=2;
    505                   md.cluster.executionpath=[pwd '/run'];
    506                   md.transient.requested_outputs={'default','MaskOceanLevelset'};
    507 
    508                   md=solveiceocean(md,'Transient','runtimename',false);
    509 
    510 %       %eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' ']);
    511 %       eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' : -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
    512 end
    513 % }}}
    514 % {{{ RunCoupledMITgcmISSMMatlabCoupling:
    515 if perform(org,'RunCoupledMITgcmISSMMatlabCoupling'),
    516 
    517         %load data:
    518         loaddata(org,'Parameters');
    519         loaddata(org,'ParameterizeIce');
    520         loaddata(org,'Bathymetry');
    521         loaddata(org,'IceSheetGeometry');
    522         endtime = round(MITgcmDeltaT * ...
    523          floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
    524 
    525         % {{{ prepare MITgcm
    526         % rename previous run directory and create new one
    527         if exist ('run.old')
    528             !\rm -rf run.old
    529         end
    530         if exist ('run')
    531             !\mv run run.old
    532         end
    533         !\mkdir run
    534         !\cp ../MITgcm/build/mitgcmuv run
    535         !\cp ../MITgcm/input/* run
    536         !\cp ../MITgcm/input/eedata_uncoupled run/eedata
    537 
    538         % initial salinity
    539         S=iniSalt*ones(Nx,Ny,Nz);
    540         writebin('run/Salt.bin',S);
    541 
    542         % initial temperature
    543         T=iniTheta*ones(Nx,Ny,Nz);
    544         writebin('run/Theta.bin',T);
    545 
    546         % initial velocity
    547         Z=zeros(Nx,Ny,Nz);
    548         writebin('run/Uvel.bin',Z);
    549         writebin('run/Vvel.bin',Z);
    550 
    551         % initial sea surface height
    552         Z=zeros(Nx,Ny);
    553         writebin('run/Etan.bin',Z);
    554 
    555         % salinity boundary conditions
    556         S=obcSalt*ones(Ny,Nz);
    557         thk=delZ*ones(Nz,1);
    558         bot=cumsum(thk);
    559         ik=find(bot<=mlDepth);
    560         S(:,ik)=mlSalt;
    561         writebin('run/OBs.bin',S);
    562 
    563         % temperature boundary conditions
    564         T=obcTheta*ones(Ny,Nz);
    565         T(:,ik)=mlTheta;
    566         writebin('run/OBt.bin',T);
    567 
    568         % zonal velocity boundary conditions
    569         U=obcUvel*ones(Ny,Nz);
    570         writebin('run/OBu.bin',U);
    571 
    572         % zero boundary conditions
    573         Z=zeros(Ny,Nz);
    574         writebin('run/zeros.bin',Z);
    575 
    576         % build parameter file data.obcs
    577         fidi=fopen('../MITgcm/input/data.obcs','r');
    578         fido=fopen('run/data.obcs','w');
    579         tline = fgetl(fidi);
    580         fprintf(fido,'%s\n',tline);
    581         while 1
    582             tline = fgetl(fidi);
    583             if ~ischar(tline), break, end
    584             %do the change here:
    585             if strcmpi(tline,' OB_Iwest = 40*1,'),
    586                 fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
    587                 continue;
    588             end
    589             if strcmpi(tline,' OB_Ieast = 40*-1,'),
    590                 fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
    591                 continue;
    592             end
    593             fprintf(fido,'%s\n',tline);
    594         end
    595         %close  files
    596         fclose(fidi);
    597         fclose(fido);
    598 
    599         %save bathymetry in MITgcm run directory
    600         writebin('run/bathymetry.bin',bathymetry);
    601         % }}}
    602 
    603         % {{{ ISSM settings:
    604 
    605         setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib')
    606         %timestepping:
    607         md.timestepping.start_time=start_time;
    608         md.timestepping.final_time=final_time;
    609         md.timestepping.time_step=time_step;
    610         md.cluster=generic('name',oshostname(),'np',2);
    611         md.results.TransientSolution.Base=md.geometry.base;
    612         md.transient.isgroundingline=1;
    613         md.transient.isthermal=0;
    614         md.groundingline.migration='SubelementMigration';
    615         md.groundingline.melt_interpolation='NoMeltOnPartiallyFloating';
    616         md.groundingline.friction_interpolation='SubelementFriction2';
    617 
    618         % }}}
    619 
    620         %start looping:
    621         results=md.results;
    622 
    623         for t=start_time:time_step:final_time
    624             disp(['Year: ' num2str(t)])
    625 
    626             %send draft from ISSM to MITgcm:
    627             draft=md.results.TransientSolution(end).Base;
    628             pos=find(md.mask.ice_levelset>0); draft(pos)=0;
    629             if t>start_time
    630                 old_draft=readbin('run/icetopo.bin',[Nx,Ny]);
    631             end
    632             writebin('run/icetopo.bin',draft);
    633 
    634             % {{{ generate MITgcm parameter file data
    635             fidi=fopen('../MITgcm/input/data','r');
    636             fido=fopen('run/data','w');
    637             tline = fgetl(fidi);
    638             fprintf(fido,'%s\n',tline);
    639             while 1
    640                 tline = fgetl(fidi);
    641                 if ~ischar(tline), break, end
    642                 %do the change here:
    643                 if strcmpi(tline,' xgOrigin = 0.0,'),
    644                     fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
    645                     continue;
    646                 end
    647                 if strcmpi(tline,' ygOrigin = -80.0,'),
    648                     fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
    649                     continue;
    650                 end
    651                 if strcmpi(tline,' delX = 20*0.25,'),
    652                     fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
    653                     continue;
    654                 end
    655                 if strcmpi(tline,' delY = 20*0.25,'),
    656                     fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
    657                     continue;
    658                 end
    659                 if strcmpi(tline,' delZ = 30*30.0,'),
    660                     fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
    661                     continue;
    662                 end
    663                 if strcmpi(tline,' endTime=2592000.,'),
    664                     fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
    665                     continue;
    666                 end
    667                 if strcmpi(tline,' deltaT=1200.0,'),
    668                     fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
    669                     continue;
    670                 end
    671                 if strcmpi(tline,' pChkptFreq=2592000.,'),
    672                     fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
    673                     continue;
    674                 end
    675                 if strcmpi(tline,' taveFreq=2592000.,'),
    676                     fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
    677                     continue;
    678                 end
    679                 if strcmpi(tline,' rhoConst=1030.,'),
    680                     fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
    681                     continue;
    682                 end
    683                 if strcmpi(tline,' rhoNil=1030.,'),
    684                     fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
    685                     continue;
    686                 end
    687                 fprintf(fido,'%s\n',tline);
    688             end
    689             %close  files
    690             fclose(fidi);
    691             fclose(fido);
    692             % }}}
    693 
    694             % {{{ generate initial MITgcm conditions
    695             ds=round(endtime/MITgcmDeltaT);
    696             if t>start_time
    697                 % Read pickup file
    698                 fnm=['run/pickup.' myint2str(ds,10) '.data'];
    699                 U=readbin(fnm,[Nx Ny Nz],1,'real*8',0);
    700                 V=readbin(fnm,[Nx Ny Nz],1,'real*8',1);
    701                 T=readbin(fnm,[Nx Ny Nz],1,'real*8',2);
    702                 S=readbin(fnm,[Nx Ny Nz],1,'real*8',3);
    703                 E=readbin(fnm,[Nx Ny],1,'real*8',8*Nz);
    704 
    705                 % find indices of locations where ice shelf retreated
    706                 h=readbin('run/hFacC.data',[Nx Ny Nz]);
    707                 msk=sum(h,3);
    708                 msk(find(msk))=1;
    709                 [iw jw]=find(msk); % horizontal indices where there is water
    710                 tmp=reshape(draft,[Nx,Ny])-old_draft;
    711                 tmp(find(tmp<0))=0;
    712                 [im jm]=find(tmp); % horizontal indices where there is melt
    713 
    714                 % Extrapolate T/S to locations where ice shelf retreated
    715                 for i=1:length(im)
    716 
    717                     % first try vertical extrapolation
    718                     in=find(h(im(i),jm(i),:));
    719                     if length(in)>0;
    720                         S(im(i),jm(i),1:min(in)  ) = S(im(i),jm(i),min(in));
    721                         T(im(i),jm(i),1:min(in)  ) = T(im(i),jm(i),min(in));
    722                         continue
    723                     end
    724 
    725                     % if not succesful, use closest neighbor horizontal extrapolation
    726                     [y c]=min((iw-im(i)).^2+(jw-jm(i)).^2);
    727                     salt=squeeze(S(iw(c),jw(c),:)); % salinity profile of closest neighbor
    728                     temp=squeeze(T(iw(c),jw(c),:)); % salinity profile of closest neighbor
    729                     in=find(h(iw(c),jw(c),:));
    730                     salt(1:min(in))=salt(min(in));
    731                     temp(1:min(in))=temp(min(in));
    732                     salt(max(in):end)=salt(max(in));
    733                     temp(max(in):end)=temp(max(in));
    734                     S(im(i),jm(i),:)=salt;
    735                     T(im(i),jm(i),:)=temp;
    736                 end
    737 
    738                 % Write initial conditions
    739                 writebin('run/Salt.bin' ,S);
    740                 writebin('run/Theta.bin',T);
    741                 writebin('run/Uvel.bin' ,U);
    742                 writebin('run/Vvel.bin' ,V);
    743                 writebin('run/Etan.bin' ,E);
    744             end
    745             % }}}
    746 
    747             % {{{ system call to run MITgcm
    748             cd run
    749             eval(['!mpiexec -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
    750             ts=round((t+time_step)*y2s/MITgcmDeltaT);
    751             eval(['!\mv STDERR.0000 STDERR_' myint2str(ts,10) '.data'])
    752             eval(['!\mv STDOUT.0000 STDOUT_' myint2str(ts,10) '.data'])
    753             eval(['!\cp hFacC.data hFacC_' myint2str(ts,10) '.data'])
    754             eval(['!\cp icetopo.bin icetopo_' myint2str(ts,10) '.data'])
    755             for fld={'S','T','U','V','Eta', ...
    756                      'SHICE_heatFluxtave','SHICE_fwFluxtave'}
    757                 eval(['!\mv ' fld{1} '.' myint2str(ds,10) '.data ' ...
    758                       fld{1} '_' myint2str(ts,10) '.data'])
    759             end
    760             cd ..
    761             % }}}
    762 
    763             %get melting rates from MITgcm
    764             %upward fresh water flux (kg/m^2/s):
    765             fnm=['run/SHICE_fwFluxtave_' myint2str(ts,10) '.data'];
    766             melting_rate=readbin(fnm,[Nx Ny]);
    767 
    768             %send averaged melting rate to ISSM
    769             %downward fresh water flux (m/y):
    770             melting_rate=-melting_rate(:)*y2s/rho_ice;
    771                  md.basalforcings.floatingice_melting_rate=melting_rate;
    772 
    773             % {{{ run ISSM and recover results
    774 
    775             md.timestepping.start_time=t;
    776             md.timestepping.final_time=t+time_step;;
    777                 md=solve(md,'Transient');
    778 
    779                 base=md.results.TransientSolution(end).Base;
    780                 thickness=md.results.TransientSolution(end).Thickness;
    781                 md.geometry.base=base;
    782                 md.geometry.thickness=thickness;
    783                 md.geometry.surface=md.geometry.base+md.geometry.thickness;
    784                 md.initialization.vx=md.results.TransientSolution(end).Vx;
    785                 md.initialization.vy=md.results.TransientSolution(end).Vy;
    786                 md.initialization.vel=md.results.TransientSolution(end).Vel;
    787                 md.initialization.pressure=md.results.TransientSolution(end).Pressure;
    788                 md.mask.ocean_levelset=md.results.TransientSolution(end).MaskOceanLevelset;
    789                 md.results.TransientSolution(end).FloatingiceMeltingRate=md.basalforcings.floatingice_melting_rate;
    790 
    791                 %save these results in the model, otherwise, they'll be wiped out
    792                 results(end+1)=md.results;
    793 
    794                 % }}}
    795 
    796 
    797         end
    798 
    799         md.results=results;
    800         savemodel(org,md);
    801 end
    802 % }}}
     374    loaddata(org,'Parameters');
     375    loaddata(org,'Bathymetry');
     376    loaddata(org,'IceSheetGeometry');
     377    loaddata(org,'ParameterizeIce');
     378    md=loadmodel(org,'RunUncoupledISSM');
     379
     380    % {{{ prepare ISSM: start from the steady-state
     381
     382    md.geometry.base=md.results.TransientSolution(end).Base;
     383    md.geometry.surface=md.results.TransientSolution(end).Surface;
     384    md.geometry.thickness=md.results.TransientSolution(end).Thickness;
     385    md.initialization.vx=md.results.TransientSolution(end).Vx;
     386    md.initialization.vy=md.results.TransientSolution(end).Vy;
     387    md.initialization.vel=md.results.TransientSolution(end).Vel;
     388    md.initialization.pressure=md.results.TransientSolution(end).Pressure;
     389    md.transient.isoceancoupling=1;
     390    md.transient.isgroundingline=0;
     391    md.masstransport.requested_outputs={'default','BasalforcingsFloatingiceMeltingRate'};
     392
     393    % }}}
     394    % {{{ prepare MITgcm
     395    % rename previous run directory and create new one
     396    if exist ('run.old')
     397        !\rm -rf run.old
     398    end
     399    if exist ('run')
     400        !\mv run run.old
     401    end
     402    !\mkdir run
     403    !\cp ../MITgcm/build/mitgcmuv run
     404    !\cp ../MITgcm/input_4003/* run
     405    !\cp ../MITgcm/input_4003/eedata_uncoupled run/eedata
     406   
     407    % {{{ Construct MITgcm binary input files
     408    namF='run/bathy_flat.bin';
     409    fid=fopen(namF,'w','b'); fwrite(fid,bathymetry,prec);fclose(fid);
     410   
     411    namF='run/shelficeTopo.Lin.bin';
     412    fid=fopen(namF,'w','b'); fwrite(fid,draft,prec);fclose(fid);
     413   
     414    var=([1:Ny]-2)/(jEnd-2);
     415    dMdt_fy=zeros(Nx,Ny);
     416    for i=1:Nx
     417        dMdt_fy(i,:)=-cos(pi*var);
     418        dMdt_fy(i,1)=0; dMdt_fy(i,j2:Ny)=0;
     419    end   
     420    %- rate of change due to ice-stream dynamics
     421    rateDyn=rho_water*0.1/3600;  sfx='r02';
     422    dMdt=rateDyn*dMdt_fy;
     423    namF=sprintf('%s.%s.%s','run/shelfice_dMdt',sfx,'bin');
     424    fid=fopen(namF,'w','b'); fwrite(fid,dMdt,prec);fclose(fid);
     425
     426    T_sfc = -1.9;
     427    T_bot = 2;
     428    del_T = (T_bot - T_sfc)/(59*delZ);
     429    tref=zeros(1,Nz);
     430    for k = 1:Nz;
     431        tref(k) = T_sfc + del_T*((k-20)*delZ);
     432        tref(k)= max(T_sfc,min(tref(k),T_bot));
     433    end
     434    namF='run/temp_obc.bin';
     435    tref=[tref; tref; tref];
     436    fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid);
     437
     438    S_sfc = 34.2;
     439    S_bot = 34.7;
     440    del_S = (S_bot - S_sfc)/(59*delZ);
     441    sref=zeros(1,Nz);
     442    for k = 1:Nz;
     443        sref(k) = S_sfc + del_S*((k-20)*delZ);
     444        sref(k)= max(S_sfc,min(sref(k),S_bot));
     445    end
     446    namF='run/salt_obc.bin';
     447    sref=[sref; sref; sref];
     448    fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid);
     449
     450    zax=[1:Nz];
     451    v1=2.5e-2;
     452    var=1+Nz-2*zax; var=var/(Nz-1);
     453    vobc=v1*var;
     454    namF='run/vVel_obc.bin';
     455    vobc=[vobc; vobc; vobc];
     456    fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid);
     457
     458    var=zeros(Nx,Ny,Nz);
     459    for i=1:Nx, for j=1:Ny
     460            var(i,j,:)=tref(1,:);
     461        end, end
     462    namF='run/temp_ini.bin';
     463    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
     464
     465    for i=1:Nx, for j=1:Ny
     466            var(i,j,:)=sref(1,:);
     467        end, end
     468    namF='run/salt_ini.bin';
     469    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
     470
     471    rhoAvr=rho_water-1.345;
     472    mIce=-rhoAvr*draft;
     473    namF='run/shelficeMass.Lin.bin';
     474    fid=fopen(namF,'w','b'); fwrite(fid,mIce,prec);fclose(fid);
     475    % }}}
     476    % }}}
     477
     478    md.timestepping.coupling_time=1/24/365;          % 1 hour in decimal years
     479    md.timestepping.time_step=1/24/365;              % 1 hour in decimal years
     480    md.timestepping.final_time=23/24/365;            % hour 23 in decimal years
     481    md.cluster.npocean=nPx*nPy;
     482    md.cluster.np=2;
     483    md.cluster.executionpath=[pwd '/run'];
     484    md.transient.requested_outputs={'default','MaskOceanLevelset'};
     485
     486    md=solveiceocean(md,'Transient','runtimename',false);
     487
     488%eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' ']);
     489%eval(['!mpiexec -np ' int2str(md.cluster.np) ' ' md.cluster.codepath '/issm_ocean.exe TransientSolution ' pwd ' ' md.miscellaneous.name ' : -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
     490end
     491% }}}
Note: See TracChangeset for help on using the changeset viewer.