confused about lmd_vmix.F

Report or discuss software problems and other woes

Moderators: arango, robertson

Post Reply
Message
Author
shiyf
Posts: 52
Joined: Mon Sep 09, 2013 11:09 pm
Location: OUC

confused about lmd_vmix.F

#1 Unread post by shiyf »

hi,
I try to modify the kpp mixing scheme through adding a coefficient to AKT and AKv, I modified some fortran program to read the bv coefficient as a forcing file,as:
nc = netcdf('E:\bvdata\bv_hrho1.nc', 'noclobber');
if isempty(nc), return, end

%% Global attributes:

%% (none)

%% Dimensions:

nc('xi_rho') = 461;
nc('eta_rho') = 237;
nc('s_rho') = 16;
nc('bv_time') = 12;

%% Variables and attributes:

nc{'s_rho'} = ncdouble('s_rho'); %% 16 elements.
nc{'s_rho'}.long_name = ncchar(''S-coordinate at RHO-points'');
nc{'s_rho'}.units = ncchar(''nondimensional'');
nc{'s_rho'}.valid_min = ncdouble(-1);
nc{'s_rho'}.valid_max = ncdouble(0);

nc{'bv_time'} = ncdouble('bv_time'); %% 12 elements.
nc{'bv_time'}.long_name = ncchar(''time since initialization'');
nc{'bv_time'}.units = ncchar(''days'');
nc{'bv_time'}.cycle_length = ncdouble(360);

nc{'lon_rho'} = ncdouble('eta_rho', 'xi_rho'); %% 109257 elements.
nc{'lon_rho'}.long_name = ncchar(''longitude of RHO-points'');
nc{'lon_rho'}.units = ncchar(''degree_east'');

nc{'lat_rho'} = ncdouble('eta_rho', 'xi_rho'); %% 109257 elements.
nc{'lat_rho'}.long_name = ncchar(''latitude of RHO-points'');
nc{'lat_rho'}.units = ncchar(''degree_north'');

nc{'BV'} = ncdouble('bv_time', 's_rho', 'eta_rho', 'xi_rho'); %% 20977344 elements.
nc{'BV'}.long_name = ncchar(''wave mixing'');
nc{'BV'}.units = ncchar(''m second-2'');
nc{'BV'}.time = ncchar(''bv_time'');

endef(nc)
close(nc)


and to be sure add this coefficient,I modify LMD_vmix.f and put out some number,as:

# ifdef WAVE_MIXING
if(i == 300..and.j==50)then
write(*,*)'k1',k,'tv1',Akv(i,j,k),Akt(i,j,k,1)
write(*,*),'bv1',bv(i,j,k)
endif

AKv(i,j,k)=AKv(i,j,k)+bv(i,j,k)*1.0_r8
AKt(i,j,k,itemp)=AKt(i,j,k,itemp)+bv(i,j,k)*1.0_r8

if(i == 300..and.j==50)then
write(*,*)'k2',k,'tv2',Akv(i,j,k),Akt(i,j,k,1)
write(*,*),'bv1',bv(i,j,k)
end
if
# endif
model is running,and I get log file show that:

STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME
C => (i,j,k) Cu Cv Cw Max Speed

0 0 00:00:00 1.142018E-07 2.095311E+04 2.095311E+04 6.117298E+17
(459,065,16) 5.075331E-04 0.000000E+00 0.000000E+00 2.215411E-01
k1 1 tv1 1.188390698940670E-003 1.188390698940670E-004
bv1 4.915707402042558E-017
k2 1 tv2 1.188390698940720E-003 1.188390698941162E-004
bv1 4.915707402042558E-017
k1 2 tv1 1.038055380879654E-003 1.038055380879654E-004
bv1 4.915707402042558E-017
k2 2 tv2 1.038055380879703E-003 1.038055380880145E-004
bv1 4.915707402042558E-017
k1 3 tv1 8.148467712667028E-004 8.148467712667028E-005
bv1 4.915707402042558E-017
k2 3 tv2 8.148467712667519E-004 8.148467712671944E-005
bv1 4.915707402042558E-017
k1 4 tv1 5.764586523517679E-004 5.764586523517680E-005
bv1 4.915707402042558E-017
k2 4 tv2 5.764586523518171E-004 5.764586523522596E-005
bv1 4.915707402042558E-017
k1 5 tv1 4.893519650554716E-004 4.893519650554716E-005
bv1 4.915707402042558E-017
k2 5 tv2 4.893519650555207E-004 4.893519650559631E-005
bv1 4.915707402042558E-017
k1 6 tv1 4.565035625271142E-004 4.565035625271142E-005
bv1 4.915707402042558E-017
k2 6 tv2 4.565035625271633E-004 4.565035625276057E-005
bv1 4.915707402042558E-017
k1 7 tv1 3.723081648698395E-004 3.723081648698395E-005
bv1 4.915707402042558E-017
k2 7 tv2 3.723081648698886E-004 3.723081648703310E-005
bv1 4.915707402042558E-017
k1 8 tv1 3.253900933909154E-004 3.253900933909154E-005
bv1 1.357094022350036E-014
k2 8 tv2 3.253900934044863E-004 3.253900935266248E-005
bv1 1.357094022350036E-014
k1 9 tv1 1.835621813462567E-004 1.835621813462567E-005
bv1 6.864456122197708E-012
k2 9 tv2 1.835621882107128E-004 1.835622499908179E-005
bv1 6.864456122197708E-012
k1 10 tv1 7.024491365679687E-005 7.024491365679687E-006
bv1 2.732198374254564E-009
k2 10 tv2 7.024764585517112E-005 7.027223564053941E-006
bv1 2.732198374254564E-009
k1 11 tv1 8.255620724961397E-005 8.255620724961397E-006
bv1 2.536174734846885E-007
k2 11 tv2 8.280982472309865E-005 8.509238198446086E-006
bv1 2.536174734846885E-007
k1 12 tv1 1.291638007646052E-004 1.291638007646052E-005
bv1 7.581283104932864E-006
k2 12 tv2 1.367450838695381E-004 2.049766318139338E-005
bv1 7.581283104932864E-006
k1 13 tv1 1.835748486414312E-004 1.835748486414312E-005
bv1 9.181330167770048E-005
k2 13 tv2 2.753881503191317E-004 1.101707865418436E-004
bv1 9.181330167770048E-005
k1 14 tv1 2.054684676595838E-004 2.054684676595838E-005
bv1 6.153457539623998E-004
k2 14 tv2 8.208142216219836E-004 6.358926007283582E-004
bv1 6.153457539623998E-004
k1 15 tv1 1.480808236878849E-004 1.480808236878849E-005
bv1 2.700192689006900E-003
k2 15 tv2 2.848273512694785E-003 2.715000771375689E-003
bv1 2.700192689006900E-003


Can be seen from the above,bv coefficient have been added to AKT and AKV;
but when i check the result, the temp,AKT,akv are no change,In a certain point, for example:

model result which no change: akt(300,50,:,1);
0
0.0011
0.0010
0.0008
0.0005
0.0005
0.0005
0.0004
0.0004
0.0002
0.0001
0.0001
0.0001
0.0004
0.0118
0.0126
0


bv coefficent:
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0007
0.0030
0.0147

model result with bv coefficient,akv(300,50,:,1):
0
0.0011
0.0009
0.0007
0.0005
0.0004
0.0004
0.0004
0.0003
0.0001
0.0001
0.0002
0.0001
0.0005
0.0093
0.0115
0

I am very confused,and struggle for along time,could you give some advice, thank you very much!

Post Reply