source: issm/oecreview/Archive/23390-24306/ISSM-23877-23878.diff

Last change on this file was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 35.2 KB

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

     
     1%Test Name: IceOceanCoupling Dan Goldberg'd setup
     2%ISSM/MITgcm coupled set-up
     3%
     4%Script control parameters
     5steps=[1 2 3 4 5 6 7 8 9 10 11 12];
     6final_time=1/365;
     7
     8%To download and recompile MITgcm from scratch:
     9rm -rf $ISSM_DIR/test/MITgcm/install
     10rm -rf $ISSM_DIR/test/MITgcm/build
     11rm $ISSM_DIR/test/MITgcm/code/SIZE.h
     12
     13%Organizer
     14mkdir Models
     15org=organizer('repository','Models/','prefix','IceOcean.','steps',steps);
     16
     17presentdirectory=pwd;
     18
     19% {{{ Parameters:
     20if perform(org,'Parameters'),
     21        Nx=3; %number of longitude cells
     22        Ny=200; %number of latitude cells
     23        Nz=90; %number of MITgcm vertical cells
     24        nPx=1; %number of MITgcm processes to use in x direction
     25        nPy=1; %number of MITgcm processes to use in y direction
     26        xgOrigin=1; %origin of longitude
     27        ygOrigin=1; %origin of latitude
     28        dLong=1; %longitude grid spacing
     29        dLat=1; %latitude grid spacing
     30        delZ=10; %thickness of vertical levels
     31        icefront_position_ratio=.75;
     32        ice_thickness=1000;
     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);
     70end
     71% }}}
     72% {{{ Bathymetry:
     73if perform(org,'Bathymetry'),
     74
     75    loaddata(org,'Parameters');
     76    %create lat,long
     77    lat=(ygOrigin+dLat/2):dLat:(ygOrigin+Ny*dLat);
     78    long=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong);
     79    [lat long]=meshgrid(lat,long);
     80
     81    longmin=min(long(:));
     82    longmax=max(long(:));
     83    latmin=min(lat(:));
     84    latmax=max(lat(:));
     85
     86    %create bedrock/bathymetry:
     87    bedrock=zeros(Nx,Ny);
     88    bedrock=hmax-deltah*tanh(pi*(2*(lat-latmin)./(latmax-latmin)-1))+ ...
     89            trough_depth*cos(2*pi*long./(longmax-longmin));
     90
     91    %save bathymetry file for MITgcm
     92    bathymetry=bedrock-sea_level;
     93    savedata(org,lat,long,bathymetry);
     94
     95end
     96% }}}
     97% {{{ IceSheetGeometry:
     98if 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);
     127end
     128% }}}
     129
     130%Configure MITgcm
     131% {{{ GetMITgcm:
     132if perform(org,'GetMITgcm'),
     133  system([pwd '/../MITgcm/get_mitgcm.sh']);
     134end
     135% }}}
     136% {{{ BuildMITgcm:
     137if perform(org,'BuildMITgcm'),
     138
     139    %load data:
     140    loaddata(org,'Parameters');
     141
     142    %specify computational grid in SIZE.h
     143    if ~exist('../MITgcm/code/SIZE.h')
     144        fidi=fopen('../MITgcm/code/SIZE.h.bak','r');
     145        fido=fopen('../MITgcm/code/SIZE.h','w');
     146        tline = fgetl(fidi);
     147        fprintf(fido,'%s\n',tline);
     148        while 1
     149            tline = fgetl(fidi);
     150            if ~ischar(tline), break, end
     151            %do the change here:
     152            if strcmpi(tline,'     &           sNx =  20,'),
     153                fprintf(fido,'%s%i%s\n','     &           sNx =  ',round(Nx/nPx),',');
     154                continue;
     155            end
     156            if strcmpi(tline,'     &           sNy =  20,'),
     157                fprintf(fido,'%s%i%s\n','     &           sNy =  ',round(Ny/nPy),',');
     158                continue;
     159            end
     160            if strcmpi(tline,'     &           nPx =   1,'),
     161                fprintf(fido,'%s%i%s\n','     &           nPx = ',nPx,',');
     162                continue;
     163            end
     164            if strcmpi(tline,'     &           nPy =   2,'),
     165                fprintf(fido,'%s%i%s\n','     &           nPy = ',nPy,',');
     166                continue;
     167            end
     168            fprintf(fido,'%s\n',tline);
     169        end
     170        %close  files
     171        fclose(fidi);
     172        fclose(fido);
     173    end
     174
     175    system(['../MITgcm/build.sh generic ' pwd '/../MITgcm']);
     176end
     177% }}}
     178addpath(recursivepath([pwd '/../MITgcm']));
     179% {{{ RunUncoupledMITgcm:
     180if perform(org,'RunUncoupledMITgcm'),
     181
     182    %load data:
     183    loaddata(org,'Parameters');
     184    loaddata(org,'Bathymetry');
     185    loaddata(org,'IceSheetGeometry');
     186         endtime = round(MITgcmDeltaT * ...
     187                 floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
     188
     189    % {{{ prepare MITgcm
     190    % rename previous run directory and create new one
     191    if exist ('run.old')
     192        !\rm -rf run.old
     193    end
     194    if exist ('run')
     195        !\mv run run.old
     196    end
     197    !\mkdir run
     198    !\cp ../MITgcm/build/mitgcmuv run
     199    !\cp ../MITgcm/input/* run
     200    !\cp ../MITgcm/input/eedata_uncoupled run/eedata
     201   
     202    %load data:
     203    loaddata(org,'Parameters');
     204
     205    % initial salinity
     206    S=iniSalt*ones(Nx,Ny,Nz);
     207    writebin('run/Salt.bin',S);
     208   
     209    % initial temperature
     210    T=iniTheta*ones(Nx,Ny,Nz);
     211    writebin('run/Theta.bin',T);
     212   
     213    % initial velocity
     214    Z=zeros(Nx,Ny,Nz);
     215    writebin('run/Uvel.bin',Z);
     216    writebin('run/Vvel.bin',Z);
     217   
     218    % initial sea surface height
     219    Z=zeros(Nx,Ny);
     220    writebin('run/Etan.bin',Z);
     221
     222    % salinity boundary conditions
     223    S=obcSalt*ones(Ny,Nz);
     224    thk=delZ*ones(Nz,1);
     225    bot=cumsum(thk);
     226    ik=find(bot<=mlDepth);
     227    S(:,ik)=mlSalt;
     228    writebin('run/OBs.bin',S);
     229
     230    % temperature boundary conditions
     231    T=obcTheta*ones(Ny,Nz);
     232    T(:,ik)=mlTheta;
     233    writebin('run/OBt.bin',T);
     234
     235    % zonal velocity boundary conditions
     236    U=obcUvel*ones(Ny,Nz);
     237    writebin('run/OBu.bin',U);
     238   
     239    % zero boundary conditions
     240    Z=zeros(Ny,Nz);
     241    writebin('run/zeros.bin',Z);
     242
     243    % build parameter file data.obcs
     244    fidi=fopen('../MITgcm/input/data.obcs','r');
     245    fido=fopen('run/data.obcs','w');
     246    tline = fgetl(fidi);
     247    fprintf(fido,'%s\n',tline);
     248    while 1
     249        tline = fgetl(fidi);
     250        if ~ischar(tline), break, end
     251        %do the change here:
     252        if strcmpi(tline,' OB_Iwest = 40*1,'),
     253            fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
     254            continue;
     255        end
     256        if strcmpi(tline,' OB_Ieast = 40*-1,'),
     257            fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
     258            continue;
     259        end
     260        fprintf(fido,'%s\n',tline);
     261    end
     262    %close  files
     263    fclose(fidi);
     264    fclose(fido);
     265
     266    %save bathymetry and bedrock in run directory
     267    writebin('run/bathymetry.bin',bathymetry);
     268    writebin('run/icetopo.bin',draft);
     269    % }}}
     270
     271    %start looping: 
     272    for t=start_time:time_step:final_time,
     273        disp(['Year: ' num2str(t)])
     274        % {{{ generate MITgcm parameter file data
     275        fidi=fopen('../MITgcm/input/data','r');
     276        fido=fopen('run/data','w');
     277        tline = fgetl(fidi);
     278        fprintf(fido,'%s\n',tline);
     279        while 1
     280            tline = fgetl(fidi);
     281            if ~ischar(tline), break, end
     282            %do the change here:
     283            if strcmpi(tline,' xgOrigin = 0.0,'),
     284                fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
     285                continue;
     286            end
     287            if strcmpi(tline,' ygOrigin = -80.0,'),
     288                fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
     289                continue;
     290            end
     291            if strcmpi(tline,' delX = 20*0.25,'),
     292                fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
     293                continue;
     294            end
     295            if strcmpi(tline,' delY = 20*0.25,'),
     296                fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
     297                continue;
     298            end
     299            if strcmpi(tline,' delZ = 30*30.0,'),
     300                fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
     301                continue;
     302            end
     303            if strcmpi(tline,' endTime=2592000.,'),
     304                fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
     305                continue;
     306            end
     307            if strcmpi(tline,' deltaT=1200.0,'),
     308                fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
     309                continue;
     310            end
     311            if strcmpi(tline,' pChkptFreq=2592000.,'),
     312                fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
     313                continue;
     314            end
     315            if strcmpi(tline,' taveFreq=2592000.,'),
     316                fprintf(fido,'%s%i%s\n',' taveFreq= ',endtime,',');
     317                continue;
     318            end
     319            if strcmpi(tline,' rhoConst=1030.,'),
     320                fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
     321                continue;
     322            end
     323            if strcmpi(tline,' rhoNil=1030.,'),
     324                fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
     325                continue;
     326            end
     327            fprintf(fido,'%s\n',tline);
     328        end
     329        %close  files
     330        fclose(fidi);
     331        fclose(fido);
     332        % }}}
     333        % {{{ generate initial MITgcm conditions
     334
     335        ds=round(endtime/MITgcmDeltaT);
     336        if t>start_time
     337            % Read pickup file
     338            fnm=['run/pickup.' myint2str(ds,10) '.data'];
     339            U=readbin(fnm,[Nx Ny Nz],1,'real*8',0);
     340            V=readbin(fnm,[Nx Ny Nz],1,'real*8',1);
     341            T=readbin(fnm,[Nx Ny Nz],1,'real*8',2);
     342            S=readbin(fnm,[Nx Ny Nz],1,'real*8',3);
     343            E=readbin(fnm,[Nx Ny],1,'real*8',8*Nz);
     344            writebin('run/Salt.bin' ,S);
     345            writebin('run/Theta.bin',T);
     346            writebin('run/Uvel.bin' ,U);
     347            writebin('run/Vvel.bin' ,V);
     348            writebin('run/Etan.bin' ,E);
     349        end
     350
     351        % }}}
     352        % {{{ system call to run MITgcm
     353        cd run
     354        eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
     355        ts=round((t+time_step)*y2s/MITgcmDeltaT);
     356        eval(['!\mv STDERR.0000 STDERR_' myint2str(ts,10) '.data'])
     357        eval(['!\mv STDOUT.0000 STDOUT_' myint2str(ts,10) '.data'])
     358        eval(['!\cp hFacC.data hFacC_' myint2str(ts,10) '.data'])
     359        eval(['!\cp icetopo.bin icetopo_' myint2str(ts,10) '.data'])
     360        for fld={'S','T','U','V','Eta', ...
     361                 'SHICE_heatFluxtave','SHICE_fwFluxtave'}
     362            eval(['!\mv ' fld{1} '.' myint2str(ds,10) '.data ' ...
     363                  fld{1} '_' myint2str(ts,10) '.data'])
     364        end
     365        cd ..
     366        % }}}
     367    end
     368end
     369% }}}
     370
     371%Configure ISSM
     372% {{{ CreateMesh:
     373if 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
     410end
     411% }}}
     412% {{{ MeshGeometry:
     413if 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);
     466end
     467% }}}
     468% {{{ ParameterizeIce:
     469if 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);
     538end
     539% }}}
     540% {{{ RunUncoupledISSM:
     541if 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);
     559end
     560% }}}
     561
     562%Run MITgcm/ISSM
     563% {{{ RunCoupledMITgcmISSM:
     564if perform(org,'RunCoupledMITgcmISSM'),
     565
     566        %load data:
     567        loaddata(org,'Parameters');
     568        loaddata(org,'ParameterizeIce');
     569        loaddata(org,'Bathymetry');
     570        loaddata(org,'IceSheetGeometry');
     571        endtime = round(MITgcmDeltaT * ...
     572         floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
     573
     574        % {{{ prepare MITgcm
     575        % rename previous run directory and create new one
     576        if exist ('run.old')
     577            !\rm -rf run.old
     578        end
     579        if exist ('run')
     580            !\mv run run.old
     581        end
     582        !\mkdir run
     583        !\cp ../MITgcm/build/mitgcmuv run
     584        !\cp ../MITgcm/input/* run
     585        !\cp ../MITgcm/input/eedata_uncoupled run/eedata
     586
     587        % initial salinity
     588        S=iniSalt*ones(Nx,Ny,Nz);
     589        writebin('run/Salt.bin',S);
     590
     591        % initial temperature
     592        T=iniTheta*ones(Nx,Ny,Nz);
     593        writebin('run/Theta.bin',T);
     594
     595        % initial velocity
     596        Z=zeros(Nx,Ny,Nz);
     597        writebin('run/Uvel.bin',Z);
     598        writebin('run/Vvel.bin',Z);
     599
     600        % initial sea surface height
     601        Z=zeros(Nx,Ny);
     602        writebin('run/Etan.bin',Z);
     603
     604        % salinity boundary conditions
     605        S=obcSalt*ones(Ny,Nz);
     606        thk=delZ*ones(Nz,1);
     607        bot=cumsum(thk);
     608        ik=find(bot<=mlDepth);
     609        S(:,ik)=mlSalt;
     610        writebin('run/OBs.bin',S);
     611
     612        % temperature boundary conditions
     613        T=obcTheta*ones(Ny,Nz);
     614        T(:,ik)=mlTheta;
     615        writebin('run/OBt.bin',T);
     616
     617        % zonal velocity boundary conditions
     618        U=obcUvel*ones(Ny,Nz);
     619        writebin('run/OBu.bin',U);
     620
     621        % zero boundary conditions
     622        Z=zeros(Ny,Nz);
     623        writebin('run/zeros.bin',Z);
     624
     625        % build parameter file data.obcs
     626        fidi=fopen('../MITgcm/input/data.obcs','r');
     627        fido=fopen('run/data.obcs','w');
     628        tline = fgetl(fidi);
     629        fprintf(fido,'%s\n',tline);
     630        while 1
     631            tline = fgetl(fidi);
     632            if ~ischar(tline), break, end
     633            %do the change here:
     634            if strcmpi(tline,' OB_Iwest = 40*1,'),
     635                fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
     636                continue;
     637            end
     638            if strcmpi(tline,' OB_Ieast = 40*-1,'),
     639                fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
     640                continue;
     641            end
     642            fprintf(fido,'%s\n',tline);
     643        end
     644        %close  files
     645        fclose(fidi);
     646        fclose(fido);
     647
     648        %save bathymetry in MITgcm run directory
     649        writebin('run/bathymetry.bin',bathymetry);
     650        % }}}
     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
     673            disp(['Year: ' num2str(t)])
     674
     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;
     678            if t>start_time
     679                old_draft=readbin('run/icetopo.bin',[Nx,Ny]);
     680            end
     681            writebin('run/icetopo.bin',draft);
     682
     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);
     688            while 1
     689                tline = fgetl(fidi);
     690                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
     708                if strcmpi(tline,' delZ = 30*30.0,'),
     709                    fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
     710                    continue;
     711                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
     728                if strcmpi(tline,' rhoConst=1030.,'),
     729                    fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
     730                    continue;
     731                end
     732                if strcmpi(tline,' rhoNil=1030.,'),
     733                    fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
     734                    continue;
     735                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
     744            ds=round(endtime/MITgcmDeltaT);
     745            if t>start_time
     746                % Read pickup file
     747                fnm=['run/pickup.' myint2str(ds,10) '.data'];
     748                U=readbin(fnm,[Nx Ny Nz],1,'real*8',0);
     749                V=readbin(fnm,[Nx Ny Nz],1,'real*8',1);
     750                T=readbin(fnm,[Nx Ny Nz],1,'real*8',2);
     751                S=readbin(fnm,[Nx Ny Nz],1,'real*8',3);
     752                E=readbin(fnm,[Nx Ny],1,'real*8',8*Nz);
     753
     754                % find indices of locations where ice shelf retreated
     755                h=readbin('run/hFacC.data',[Nx Ny Nz]);
     756                msk=sum(h,3);
     757                msk(find(msk))=1;
     758                [iw jw]=find(msk); % horizontal indices where there is water
     759                tmp=reshape(draft,[Nx,Ny])-old_draft;
     760                tmp(find(tmp<0))=0;
     761                [im jm]=find(tmp); % horizontal indices where there is melt
     762               
     763                % Extrapolate T/S to locations where ice shelf retreated
     764                for i=1:length(im)
     765
     766                    % first try vertical extrapolation
     767                    in=find(h(im(i),jm(i),:));
     768                    if length(in)>0;
     769                        S(im(i),jm(i),1:min(in)  ) = S(im(i),jm(i),min(in));
     770                        T(im(i),jm(i),1:min(in)  ) = T(im(i),jm(i),min(in));
     771                        continue
     772                    end
     773
     774                    % if not succesful, use closest neighbor horizontal extrapolation
     775                    [y c]=min((iw-im(i)).^2+(jw-jm(i)).^2);
     776                    salt=squeeze(S(iw(c),jw(c),:)); % salinity profile of closest neighbor
     777                    temp=squeeze(T(iw(c),jw(c),:)); % salinity profile of closest neighbor
     778                    in=find(h(iw(c),jw(c),:));
     779                    salt(1:min(in))=salt(min(in));
     780                    temp(1:min(in))=temp(min(in));
     781                    salt(max(in):end)=salt(max(in));
     782                    temp(max(in):end)=temp(max(in));
     783                    S(im(i),jm(i),:)=salt;
     784                    T(im(i),jm(i),:)=temp;
     785                end
     786
     787                % Write initial conditions
     788                writebin('run/Salt.bin' ,S);
     789                writebin('run/Theta.bin',T);
     790                writebin('run/Uvel.bin' ,U);
     791                writebin('run/Vvel.bin' ,V);
     792                writebin('run/Etan.bin' ,E);
     793            end
     794            % }}}
     795
     796            % {{{ system call to run MITgcm
     797            cd run
     798            eval(['!mpiexec -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
     799            ts=round((t+time_step)*y2s/MITgcmDeltaT);
     800            eval(['!\mv STDERR.0000 STDERR_' myint2str(ts,10) '.data'])
     801            eval(['!\mv STDOUT.0000 STDOUT_' myint2str(ts,10) '.data'])
     802            eval(['!\cp hFacC.data hFacC_' myint2str(ts,10) '.data'])
     803            eval(['!\cp icetopo.bin icetopo_' myint2str(ts,10) '.data'])
     804            for fld={'S','T','U','V','Eta', ...
     805                     'SHICE_heatFluxtave','SHICE_fwFluxtave'}
     806                eval(['!\mv ' fld{1} '.' myint2str(ds,10) '.data ' ...
     807                      fld{1} '_' myint2str(ts,10) '.data'])
     808            end
     809            cd ..
     810            % }}}
     811
     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);
     850end
     851% }}}
     852% {{{ RunCoupledMITgcmISSM2:
     853if 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
     863        % rename previous run directory and create new one
     864        if exist ('run.old')
     865            !\rm -rf run.old
     866        end
     867        if exist ('run')
     868            !\mv run run.old
     869        end
     870        !\mkdir run
     871        !\cp ../MITgcm/build/mitgcmuv run
     872        !\cp ../MITgcm/input/* run
     873
     874        % initial salinity
     875        S=iniSalt*ones(Nx,Ny,Nz);
     876        writebin('run/Salt.bin',S);
     877
     878        % initial temperature
     879        T=iniTheta*ones(Nx,Ny,Nz);
     880        writebin('run/Theta.bin',T);
     881
     882        % initial velocity
     883        Z=zeros(Nx,Ny,Nz);
     884        writebin('run/Uvel.bin',Z);
     885        writebin('run/Vvel.bin',Z);
     886
     887        % initial sea surface height
     888        Z=zeros(Nx,Ny);
     889        writebin('run/Etan.bin',Z);
     890
     891        % salinity boundary conditions
     892        S=obcSalt*ones(Ny,Nz);
     893        thk=delZ*ones(Nz,1);
     894        bot=cumsum(thk);
     895        ik=find(bot<=mlDepth);
     896        S(:,ik)=mlSalt;
     897        writebin('run/OBs.bin',S);
     898
     899        % temperature boundary conditions
     900        T=obcTheta*ones(Ny,Nz);
     901        T(:,ik)=mlTheta;
     902        writebin('run/OBt.bin',T);
     903
     904        % zonal velocity boundary conditions
     905        U=obcUvel*ones(Ny,Nz);
     906        writebin('run/OBu.bin',U);
     907
     908        % zero boundary conditions
     909        Z=zeros(Ny,Nz);
     910        writebin('run/zeros.bin',Z);
     911
     912        % build parameter file data.obcs
     913        fidi=fopen('../MITgcm/input/data.obcs','r');
     914        fido=fopen('run/data.obcs','w');
     915        tline = fgetl(fidi);
     916        fprintf(fido,'%s\n',tline);
     917        while 1
     918            tline = fgetl(fidi);
     919            if ~ischar(tline), break, end
     920            %do the change here:
     921            if strcmpi(tline,' OB_Iwest = 40*1,'),
     922                fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
     923                continue;
     924            end
     925            if strcmpi(tline,' OB_Ieast = 40*-1,'),
     926                fprintf(fido,'%s%i%s\n',' OB_Ieast = ',Ny,'*-1,');
     927                continue;
     928            end
     929            fprintf(fido,'%s\n',tline);
     930        end
     931        %close  files
     932        fclose(fidi);
     933        fclose(fido);
     934
     935        %save bathymetry and bedrock in run directory
     936        writebin('run/bathymetry.bin',bathymetry);
     937        writebin('run/icetopo.bin',draft);
     938        % }}}
     939        % {{{ generate MITgcm parameter file data
     940        fidi=fopen('../MITgcm/input/data','r');
     941        fido=fopen('run/data','w');
     942        tline = fgetl(fidi);
     943        fprintf(fido,'%s\n',tline);
     944        while 1
     945            tline = fgetl(fidi);
     946            if ~ischar(tline), break, end
     947            %do the change here:
     948            if strcmpi(tline,' xgOrigin = 0.0,'),
     949                fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
     950                continue;
     951            end
     952            if strcmpi(tline,' ygOrigin = -80.0,'),
     953                fprintf(fido,'%s%i%s\n',' ygOrigin = ',ygOrigin,',');
     954                continue;
     955            end
     956            if strcmpi(tline,' delX = 20*0.25,'),
     957                fprintf(fido,'%s%i*%g%s\n',' delX = ',Nx,dLong,',');
     958                continue;
     959            end
     960            if strcmpi(tline,' delY = 20*0.25,'),
     961                fprintf(fido,'%s%i*%g%s\n',' delY = ',Ny,dLat,',');
     962                continue;
     963            end
     964            if strcmpi(tline,' delZ = 30*30.0,'),
     965                fprintf(fido,'%s%i*%g%s\n',' delZ = ',Nz,delZ,',');
     966                continue;
     967            end
     968            if strcmpi(tline,' endTime=2592000.,'),
     969                fprintf(fido,'%s%i%s\n',' endTime= ',endtime,',');
     970                continue;
     971            end
     972            if strcmpi(tline,' deltaT=1200.0,'),
     973                fprintf(fido,'%s%i%s\n',' deltaT= ',MITgcmDeltaT,',');
     974                continue;
     975            end
     976            if strcmpi(tline,' pChkptFreq=2592000.,'),
     977                fprintf(fido,'%s%i%s\n',' pChkptFreq= ',endtime,',');
     978                continue;
     979            end
     980            if strcmpi(tline,' taveFreq=2592000.,'),
     981                fprintf(fido,'%s%i%s\n',' taveFreq= ',outputtime,',');
     982                continue;
     983            end
     984            if strcmpi(tline,' rhoConst=1030.,'),
     985                fprintf(fido,'%s%i%s\n',' rhoConst= ',rho_water,',');
     986                continue;
     987            end
     988            if strcmpi(tline,' rhoNil=1030.,'),
     989                fprintf(fido,'%s%i%s\n',' rhoNil= ',rho_water,',');
     990                continue;
     991            end
     992            fprintf(fido,'%s\n',tline);
     993        end
     994        %close  files
     995        fclose(fidi);
     996        fclose(fido);
     997        % }}}
     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;
     1005   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']);
     1016end
     1017% }}}
     1018
     1019%Fields and tolerances to track changes
     1020fnm=['run/SHICE_fwFluxtave.0000000006.data'];
     1021melting_rate_1=readbin(fnm,[Nx Ny]);
     1022fnm=['run/SHICE_fwFluxtave.0000000012.data'];
     1023melting_rate_2=readbin(fnm,[Nx Ny]);
     1024fnm=['run/SHICE_fwFluxtave.0000000018.data'];
     1025melting_rate_3=readbin(fnm,[Nx Ny]);
     1026fnm=['run/SHICE_fwFluxtave.0000000024.data'];
     1027melting_rate_4=readbin(fnm,[Nx Ny]);
     1028field_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'};
     1031field_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 };
     1034field_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        };
Note: See TracBrowser for help on using the repository browser.