Changeset 26664


Ignore:
Timestamp:
11/24/21 13:34:49 (3 years ago)
Author:
dmenemen
Message:

test4003 now runs MITgcm/verification/shelfice_2d_remesh
except that model domain is 3x200 instead of 1x200 and
initial conditions are from time step 0 so no pickup
files are needed. Coupling with ISSM is work in progress.

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/test/MITgcm/code_4003/CPP_OPTIONS.h

    r26544 r26664  
    5454#define INCLUDE_PHIHYD_CALCULATION_CODE
    5555
    56 C o Include/exclude sound speed calculation code
    57 C o (Note that this is a diagnostic from Del Grasso algorithm, not derived from EOS)
    58 #undef INCLUDE_SOUNDSPEED_CALC_CODE
    59 
    6056C-- Vertical mixing code options:
    6157
    … …  
    7874#undef EXCLUDE_PCELL_MIX_CODE
    7975
    80 C o Exclude/allow to use isotropic 3-D Smagorinsky viscosity as diffusivity
    81 C   for tracers (after scaling by constant Prandtl number)
    82 #undef ALLOW_SMAG_3D_DIFFUSIVITY
    83 
    8476C-- Time-stepping code options:
    8577
    … …  
    9385#undef ALLOW_ADAMSBASHFORTH_3
    9486
    95 C o Include/exclude Quasi-Hydrostatic Stagger Time-step AdamsBashforth code
    96 #undef ALLOW_QHYD_STAGGER_TS
    97 
    9887C-- Model formulation options:
    9988
    … …  
    10594C   this implies that grid-cell thickness (hFactors) varies with time
    10695#define NONLIN_FRSURF
    107 C o Disable code for rStar coordinate and/or code for Sigma coordinate
    108 c#define DISABLE_RSTAR_CODE
    109 c#define DISABLE_SIGMA_CODE
    11096
    11197C o Include/exclude nonHydrostatic code
    … …  
    117103C-- Algorithm options:
    118104
    119 C o Include/exclude code for Non Self-Adjoint (NSA) conjugate-gradient solver
     105C o Use Non Self-Adjoint (NSA) conjugate-gradient solver
    120106#undef ALLOW_CG2D_NSA
    121107
    … …  
    131117C-- Retired code options:
    132118
    133 C o ALLOW isotropic scaling of harmonic and bi-harmonic terms when
    134 C   using an locally isotropic spherical grid with (dlambda) x (dphi*cos(phi))
    135 C *only for use on a lat-lon grid*
    136 C   Setting this flag here affects both momentum and tracer equation unless
    137 C   it is set/unset again in other header fields (e.g., GAD_OPTIONS.h).
    138 C   The definition of the flag is commented to avoid interference with
    139 C   such other header files.
    140 C   The preferred method is specifying a value for viscAhGrid or viscA4Grid
    141 C   in data which is then automatically scaled by the grid size;
    142 C   the old method of specifying viscAh/viscA4 and this flag is provided
    143 C   for completeness only (and for use with the adjoint).
    144 c#define ISOTROPIC_COS_SCALING
    145 
    146 C o This flag selects the form of COSINE(lat) scaling of bi-harmonic term.
    147 C *only for use on a lat-lon grid*
    148 C   Has no effect if ISOTROPIC_COS_SCALING is undefined.
    149 C   Has no effect on vector invariant momentum equations.
    150 C   Setting this flag here affects both momentum and tracer equation unless
    151 C   it is set/unset again in other header fields (e.g., GAD_OPTIONS.h).
    152 C   The definition of the flag is commented to avoid interference with
    153 C   such other header files.
    154 c#define COSINEMETH_III
    155 
    156 C o Use "OLD" UV discretisation near boundaries (*not* recommended)
    157 C   Note - only works with pkg/mom_fluxform and "no_slip_sides=.FALSE."
    158 C          because the old code did not have no-slip BCs
    159 #undef OLD_ADV_BCS
    160 
    161119C o Use LONG.bin, LATG.bin, etc., initialization for ini_curviliear_grid.F
    162120C   Default is to use "new" grid files (OLD_GRID_IO undef) but OLD_GRID_IO
    163121C   is still useful with, e.g., single-domain curvilinear configurations.
    164122#undef OLD_GRID_IO
    165 
    166 C o Use old EXTERNAL_FORCING_U,V,T,S subroutines (for backward compatibility)
    167 #undef USE_OLD_EXTERNAL_FORCING
    168123
    169124C-- Other option files:

  • issm/trunk-jpl/test/MITgcm/code_4003/DIAGNOSTICS_SIZE.h

    r26544 r26664  
    2323      PARAMETER( diagSt_size = 10*Nr )
    2424
    25 
    2625CEH3 ;;; Local Variables: ***
    2726CEH3 ;;; mode:fortran ***

  • issm/trunk-jpl/test/MITgcm/code_4003/SIZE.h

    r26545 r26664  
    4343      INTEGER Nr
    4444      PARAMETER (
    45      &           sNx =   1,
    46      &           sNy =  50,
     45     &           sNx =   3,
     46     &           sNy =  25,
    4747     &           OLx =   3,
    4848     &           OLy =   3,
    … …  
    5050     &           nSy =   1,
    5151     &           nPx =   1,
    52      &           nPy =   4,
     52     &           nPy =   8,
    5353     &           Nx  = sNx*nSx*nPx,
    5454     &           Ny  = sNy*nSy*nPy,

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

    r26545 r26664  
    6262#Time stepping parameters
    6363 &PARM03
    64 #nIter0=2880,
    65 #nTimeSteps=288,
    66 #startTime=0.,
     64 startTime=0.,
     65 nTimeSteps=192,
    6766#endTime=2592000.,
    68  deltaT=300.0,
     67 deltaT=450.0,
    6968 forcing_In_AB  = .FALSE.,
    7069 abEps=0.1,
    … …  
    7271 chkptFreq = 2592000.,
    7372 dumpFreq =  86400.,
    74  monitorFreq= 7200.,
    75  monitorSelect=1,
    76  monitorFreq=1800.,
    77  dumpFreq = 10800.,
    78 #- for testing purpose:
    79  nIter0=2898,
    80  nTimeSteps=20,
    81 #dumpFreq = 1200.,
    8273 monitorFreq=1.,
    8374 &
    … …  
    8879 usingSphericalPolarGrid=.TRUE.,
    8980 delR=90*10.,
    90  delX=1*.125,
     81 delX=3*.125,
    9182 delY=200*.0078125,
    9283 xgOrigin = 0.,

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

    r26614 r26664  
    55steps=1:12;
    66steps=1:6;
    7 steps=7;
    87final_time=1/365;
    98
    109%To download and recompile MITgcm from scratch:
    11 %!rm -rf ${ISSM_DIR}/test/MITgcm/install
    12 %!rm -rf ${ISSM_DIR}/test/MITgcm/build/*
    13 %!rm -rf Models
     10!rm -rf ${ISSM_DIR}/test/MITgcm/install
     11!rm -rf ${ISSM_DIR}/test/MITgcm/build/*
     12!rm -rf Models
    1413
    1514%Organizer
    … …  
    2120% {{{ Parameters:
    2221if perform(org,'Parameters'),
    23         Nx=1; %number of longitude cells
     22        Nx=3; %number of longitude cells
    2423        Ny=200; %number of latitude cells
    2524        Nz=90; %number of MITgcm vertical cells
    2625        nPx=1; %number of MITgcm processes to use in x direction
    27         nPy=4; %number of MITgcm processes to use in y direction
    28         xgOrigin=1; %origin of longitude
    29         ygOrigin=1; %origin of latitude
    30         dLong=1; %longitude grid spacing
    31         dLat=1; %latitude grid spacing
     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
    3231        delZ=10; %thickness of vertical levels
    3332        icefront_position_ratio=.75;
    … …  
    4746        obcUvel  = -0.1; % open boundary velocity (m/s)
    4847
    49         MITgcmDeltaT=600; % MITgcm time step in seconds
     48        MITgcmDeltaT=450; % MITgcm time step in seconds
    5049        y2s=31536000; % year to seconds conversion, i.e., seconds per year
    5150
    … …  
    5554        async_step_MITgcm_multiplier=1; % used to reduce run time for MItgcm
    5655
    57         % bedrock/bathymetry
     56        % bathymetry
    5857        hmax=1000;
    5958        trough_depth=200;
    … …  
    7675
    7776    loaddata(org,'Parameters');
    78     %create lat,long
     77
     78    v0 = 2e3;
     79    h0 = 800;
     80
     81    hfacMin = 0.2;
     82
     83    eos = 'jmd95z';
     84    prec = 'real*8';
     85
     86    %create lat,lon
     87    latg = ygOrigin+[0:Ny-1]*dLat;
     88    latc = latg+dLat/2;
     89    long = xgOrigin+[0:Nx-1]*dLong;
     90    lonc = long+dLong/2;
     91    zC=-delZ*([1:Nz]-0.5);
     92    zF=-delZ*[0:Nz];
     93           
    7994    lat=(ygOrigin+dLat/2):dLat:(ygOrigin+Ny*dLat);
    80     long=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong);
    81     [lat long]=meshgrid(lat,long);
    82 
    83     longmin=min(long(:));
    84     longmax=max(long(:));
     95    lon=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong);
     96    [lat lon]=meshgrid(lat,lon);
     97
     98    longmin=min(lon(:));
     99    longmax=max(lon(:));
    85100    latmin=min(lat(:));
    86101    latmax=max(lat(:));
    87102
    88     %create bedrock/bathymetry:
    89     bedrock=zeros(Nx,Ny);
    90     bedrock=hmax-deltah*tanh(pi*(2*(lat-latmin)./(latmax-latmin)-1))+ ...
    91             trough_depth*cos(2*pi*long./(longmax-longmin));
     103% Gravity
     104gravity= 9.81;
     105rhoConst= 1030;
     106
     107% Nominal depth of model (meters)
     108H = -900;               %water depth in the ice shelf cavity
     109Hmin = -600;            % deepest point of cavern               
     110Hmax = -300;            % shallowest point of cavern
     111jEnd = Ny*3/4;           % where ice-shelf ends
     112dHdy = (Hmax-Hmin)/dLat/(jEnd-2); %Slope of ice shelf
     113
     114    %create bathymetry:
     115    bathymetry = ones(Nx,Ny)*H;
     116    bathymetry(:,1) = 0;
    92117
    93118    %save bathymetry file for MITgcm
    94     bathymetry=bedrock-sea_level;
    95     savedata(org,lat,long,bathymetry);
     119    savedata(org,lat,lon,bathymetry);
     120
     121j2=jEnd+1;
     122hIce=bathymetry;
     123for i=1:Nx
     124    hIce(i,:)=Hmin+dHdy*[-1:Ny-2]*dLat;
     125    hIce(i,1)=0; hIce(i,j2:Ny)=0;
     126end
     127
     128var=([1:Ny]-2)/(jEnd-2);
     129dMdt_fy=bathymetry;
     130for i=1:Nx
     131    dMdt_fy(i,:)=-cos(pi*var);
     132    dMdt_fy(i,1)=0; dMdt_fy(i,j2:Ny)=0;
     133end
     134
     135regMsk=ones(Nx,Ny);
     136regMsk(:,1)=0; regMsk(:,j2:Ny)=2;
     137
     138%- rate of change due to ice-stream dynamics
     139rateDyn=rhoConst*0.1/3600;  sfx='r02';
     140
     141dMdt=rateDyn*dMdt_fy;
     142
     143dz = delZ*ones(1,Nz);
     144zgp1 = [0,cumsum(dz)];
     145zc = .5*(zgp1(1:end-1)+zgp1(2:end));
     146zg = zgp1(1:end-1);
     147dz = diff(zgp1);
     148
     149T_sfc = -1.9;
     150T_bot = 2;
     151del_T = (T_bot - T_sfc)/(59*delZ);
     152tref=0*dz;
     153for k = 1:Nz;
     154    tref(k) = T_sfc + del_T*((k-20)*delZ);
     155    tref(k)= max(T_sfc,min(tref(k),T_bot));
     156end
     157
     158S_sfc = 34.2;
     159S_bot = 34.7;
     160del_S = (S_bot - S_sfc)/(59*delZ);
     161sref=0*dz;
     162for k = 1:Nz;
     163    sref(k) = S_sfc + del_S*((k-20)*delZ);
     164    sref(k)= max(S_sfc,min(sref(k),S_bot));
     165end
     166pEOS=-rhoConst*gravity*zC; % in Pa
     167pEOS=pEOS*1.e-4; % in dBar
     168rhoAn=densjmd95(sref,tref,pEOS);
     169rhoAn=rhoAn-rhoConst;
     170
     171pF=-rhoConst*gravity*zF*1.e-4; % in dBar
     172rhoUp=densjmd95(sref,tref,pF(2:end));
     173rhoDw=densjmd95(sref,tref,pF(1:Nz));
     174dRho=rhoUp(1:Nz-1)-rhoDw(2:Nz);
     175NSq=-gravity*dRho/delZ/rhoConst;
     176
     177mnV=min(NSq); MxV=max(NSq); Avr=mean(NSq);
     178
     179zax=[1:Nz];
     180
     181v1=2.5e-2;
     182var=1+Nz-2*zax; var=var/(Nz-1);
     183vobc=v1*var;
     184
     185rhoAvr=rhoConst-1.345;
     186mIce=-rhoAvr*hIce;
    96187
    97188end
    … …  
    153244    !\mkdir run
    154245    !\cp ../MITgcm/build/mitgcmuv run
    155     !\cp ../MITgcm/install/verification/shelfice_2d_remesh/input/* run
    156246    !\cp ../MITgcm/input_4003/* run
    157     !\cp ../MITgcm/input_4003/data_uncoupled run/data
    158247    !\cp ../MITgcm/input_4003/eedata_uncoupled run/eedata
    159248    cd run
     249   
     250    namF='bathy_flat.bin';
     251    fid=fopen(namF,'w','b'); fwrite(fid,bathymetry,prec);fclose(fid);
     252   
     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');
     260    fid=fopen(namF,'w','b'); fwrite(fid,dMdt,prec);fclose(fid);
     261
     262    namF='temp_obc.bin';
     263    tref=[tref; tref; tref];
     264    fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid);
     265
     266    namF='salt_obc.bin';
     267    sref=[sref; sref; sref];
     268    fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid);
     269
     270    namF='vVel_obc.bin';
     271    vobc=[vobc; vobc; vobc];
     272    fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid);
     273
     274    var=zeros(Nx,Ny,Nz);
     275    for i=1:Nx, for j=1:Ny
     276            var(i,j,:)=tref(1,:);
     277        end, end
     278    namF='temp_ini.bin';
     279    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
     280
     281    for i=1:Nx, for j=1:Ny
     282            var(i,j,:)=sref(1,:);
     283        end, end
     284    namF='salt_ini.bin';
     285    fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid);
     286
     287    namF='shelficeMass.Lin.bin';
     288    fid=fopen(namF,'w','b'); fwrite(fid,mIce,prec);fclose(fid);
     289
    160290    eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
    161         cd ..
     291    cd ..
    162292end
    163293% }}}
    … …  
    671801end
    672802% }}}
    673 
    674 %Fields and tolerances to track changes
    675 %fnm=['run/SHICE_fwFlux.0000011400.data'];
    676 %melting_rate_1=readbin(fnm,[Nx Ny]);
    677 %fnm=['run/SHICE_fwFlux.0000011450.data'];
    678 %melting_rate_2=readbin(fnm,[Nx Ny]);
    679 %fnm=['run/SHICE_fwFlux.0000011500.data'];
    680 %melting_rate_3=readbin(fnm,[Nx Ny]);
    681 %fnm=['run/SHICE_fwFlux.0000011550.data'];
    682 %melting_rate_4=readbin(fnm,[Nx Ny]);
    683 %field_names     ={'Base1','Melting1','Vx2','Vy2','Thickness2','Base2','MaskOceanLevelset2','FloatingiceMeltingRate2',...
    684 %       'Melting2','Vx3','Vy3','Thickness3','Base3','MaskOceanLevelset3','FloatingiceMeltingRate3',...
    685 %       'Melting3','Vx4','Vy4','Thickness4','Base4','MaskOceanLevelset4','FloatingiceMeltingRate4','Melting4'};
    686 %field_tolerances={2e-13,1e-13,7e-13,7e-13,1e-13,1e-13,1e-13,1e-12,...
    687 %       1e-13, 1e-12, 1e-12, 1e-13, 1e-13, 1e-13, 1e-12,...
    688 %       1e-13, 1e-12, 1e-12, 1e-13, 1e-13, 1e-13, 1e-12, 1e-13 };
    689 %field_values={...
    690 %       (md.results.TransientSolution(1).Base),...
    691 %       (melting_rate_1(:)),...
    692 %       (md.results.TransientSolution(2).Vx),...
    693 %       (md.results.TransientSolution(2).Vy),...
    694 %       (md.results.TransientSolution(2).Thickness),...
    695 %       (md.results.TransientSolution(2).Base),...
    696 %       (md.results.TransientSolution(2).MaskOceanLevelset),...
    697 %       (md.results.TransientSolution(2).BasalforcingsFloatingiceMeltingRate),...
    698 %       (melting_rate_2(:)),...
    699 %       (md.results.TransientSolution(3).Vx),...
    700 %       (md.results.TransientSolution(3).Vy),...
    701 %       (md.results.TransientSolution(3).Thickness),...
    702 %       (md.results.TransientSolution(3).Base),...
    703 %       (md.results.TransientSolution(3).MaskOceanLevelset),...
    704 %       (md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),...
    705 %       (melting_rate_3(:)),...
    706 %       (md.results.TransientSolution(4).Vx),...
    707 %       (md.results.TransientSolution(4).Vy),...
    708 %       (md.results.TransientSolution(4).Thickness),...
    709 %       (md.results.TransientSolution(4).Base),...
    710 %       (md.results.TransientSolution(4).MaskOceanLevelset),...
    711 %       (md.results.TransientSolution(4).BasalforcingsFloatingiceMeltingRate),...
    712 %       (melting_rate_4(:)),...
    713 %       };
Note: See TracChangeset for help on using the changeset viewer.