Skip to content

Commit c0ae6c2

Browse files
jmaerzmatsbnTomasTorsvikmilicakJorgSchwinger
authored
Merge latest master into feature-hamocc_beyond-CMIP6 branch (NorESMhub#232)
* Dynamic mapping of pore water tracers to ocean tracers (NorESMhub#192) * Initial restructuring of sediment-related tracer declaration and initialization * Introducing mapping function * Remove unncessary comments * Fixed diagnostics bug and updated index naming * Added initial support for NUOPC driver. * Lon-lat variable sediment porosity (NorESMhub#189) Introducing a static 3D sediment porosity field that can be optionally read in with effects on molecular pore water diffusion and shifting. * Added wave forcing fields. * Renamed folder for MCT driver. * Moved MCT specific file from drivers/cpl_share/ to drivers/mct/. * Rename drivers/mct/mod_swtfrz.F to drivers/mct/mod_swtfrz.F90. * Rewrite to drivers/mct/mod_swtfrz.F90 to free format Fortran. * Remove redundant definition of kOBL. * Redefine kOBL, cast as integer * Fixing variable sediment porosity - field initialization in case of `sedbypass=true` (NorESMhub#198) * Removing bodensed - Initialization of sediment parameters and fields now in mo_sedmnt * This is the first commit of MKS units. All variables in the subroutines have been converted to MKS [meter, kg, seconds] instead of CGS [cm, gram, seconds] with necessary coefficients. The default option which is CGS reproduce old results. The new option MKS cannot reproduce because of machine precision. * Hamocc hybrid coord2 (NorESMhub#179) Make the surface mixed layer depth fractional index `hOBL` available for use in iHAMOCC, and adjust the internal iHAMOCC index `kmle` according to `hOBL`. Default value `kmle = 2` is retained for consistency with isopycnic coordinates. * BLOM CIME cpp updates to run in NorESM * bug fixes for the CGS MKS conversion * cesm thermal forcing bug fixes for reproducibility * BLOM MKS update to export winds into the CESM using proper units. * input values in ocn_in case is updated for mks setup * default cgsmks value changed * Initialize some variables in the k-epsilon model. * Fix porosity read (NorESMhub#201) * Fixing the reading of variable porosity input field in preparation for the NorESM 2.0.6 release Cherry-picked from private Ncycleprivate branch 0d56930e2fdd62caba964d375b57304942568926 * Provide number of layers (3rd dim) via ks and not hard-coded * minor clean-up * Correct unit of diagnostic variable dp_trc. * Made conservation and checksum diagnostics selectable by namelist options (default off). * pCO2, Piston velocity and solubility output (NorESMhub#202) * add pCO2m (moist), CO2 piston velocity and solubility output - caution: kwco2 piston velocity now really holds only piston velocity (and not times solubility) * Bugfix pnetcdf (NorESMhub#208) * Add variables used by PNETCDF to explicit use staements. * Move implicit none statments * update explicit use statement for pnetcdf * fixed units and renamed calcium burial to CaCO3 burial (NorESMhub#212) Fixed sediment clay units. * - Made the "fuk95" configuration work with MKS units. - Removed "CGS" CPP flag. - Changed some unit conversion factors from variables to parameters. - Introduced rho0 (= 1/alpha0) parameter. - Updated copyright statements. * Correct unit conversion of mixed layer depth to pressure. * Updated NorESM coupling scripts for the use of MKS units. * Fixed check of unit system when building as NorESM component. * Add option for surface pH output (NorESMhub#221) * Remove unused parameters in wrt* subroutine calls in ncout_hamocc.F90 * Import get_bgc_namelist only in subroutine where it is needed. (NorESMhub#225) * fix missing ' (NorESMhub#228) Fixing a missing ' that only showed up when using `cisonew` --------- Co-authored-by: Mats Bentsen <mben@norceresearch.no> Co-authored-by: Tomas Torsvik <tomas.torsvik.work@gmail.com> Co-authored-by: Mehmet Ilicak <milicak@itu.edu.tr> Co-authored-by: Tomas Torsvik <43031053+TomasTorsvik@users.noreply.github.com> Co-authored-by: Jörg Schwinger <jorg.schwinger@norceresearch.no>
1 parent 26e65d6 commit c0ae6c2

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

43 files changed

+828
-557
lines changed

‎ben02/mod_ben02.F

+15-6
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
! ------------------------------------------------------------------------------
2-
! Copyright (C) 2002-2021 Mats Bentsen, Mehmet Ilicak
2+
! Copyright (C) 2002-2022 Mats Bentsen, Mehmet Ilicak
33
!
44
! This file is part of BLOM.
55
!
@@ -26,7 +26,7 @@ module mod_ben02
2626
c
2727
use mod_types, only: i2, r4
2828
use mod_config, only: expcnf
29-
use mod_constants, only: t0deg, spval
29+
use mod_constants, only: t0deg, spval, L_mks2cgs
3030
use mod_calendar, only: date_offset, calendar_noerr,
3131
. calendar_errstr
3232
use mod_time, only: date, calendar, nday_in_year, nday_of_year,
@@ -183,10 +183,17 @@ module mod_ben02
183183
. atm_cswa_era ! short-wave radiation adjustment factor
184184
! (NCEP)
185185
c
186+
#ifdef MKS
187+
data atm_ice_csmt_ncep,atm_rnf_csmt_ncep /2.e10,1.e9/,
188+
. atm_crnf_ncep,atm_cswa_ncep /0.82073,0.88340/,
189+
. atm_ice_csmt_era,atm_rnf_csmt_era /0.0,1.e9/,
190+
. atm_crnf_era,atm_cswa_era /0.7234,0.9721/
191+
#else
186192
data atm_ice_csmt_ncep,atm_rnf_csmt_ncep /2.e14,1.e13/,
187193
. atm_crnf_ncep,atm_cswa_ncep /0.82073,0.88340/,
188194
. atm_ice_csmt_era,atm_rnf_csmt_era /0.0,1.e13/,
189195
. atm_crnf_era,atm_cswa_era /0.7234,0.9721/
196+
#endif
190197
c
191198
real ::
192199
. zu, ! measurement height of wind [m]
@@ -2090,11 +2097,13 @@ subroutine inifrc_ben02clim
20902097
integer, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,12) :: smtmsk
20912098
real dx2,dy2,prc_sum,eva_sum,rnf_sum,swa_sum,lwa_sum,lht_sum,
20922099
. sht_sum,fwf_fac,dangle,garea,le,albedo,fac,swa_ave,lwa_ave,
2093-
. lht_ave,sht_ave,crnf,cswa
2100+
. lht_ave,sht_ave,crnf,cswa,A_cgs2mks
20942101
real*4 rw4
20952102
integer i,j,k,l,il,jl
20962103
integer*2 rn2,ri2,rj2
20972104
c
2105+
A_cgs2mks=1./(L_mks2cgs**2)
2106+
c
20982107
c --- Allocate memory for additional monthly forcing fields.
20992108
allocate(taud (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,12),
21002109
. tauxd (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,12),
@@ -2775,7 +2784,7 @@ subroutine inifrc_ben02clim
27752784
do k=1,12
27762785
do l=1,isp(j)
27772786
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
2778-
garea=scp2(i,j)*1.e-4 ! [m^2]
2787+
garea=scp2(i,j)*A_cgs2mks ! [m^2]
27792788
c
27802789
c --- ----- freshwater fluxes [m/s]
27812790
util1(i,j)=util1(i,j)+precip(i,j,k)*fwf_fac*garea
@@ -2819,7 +2828,7 @@ subroutine inifrc_ben02clim
28192828
do j=1,jj
28202829
do l=1,isp(j)
28212830
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
2822-
garea=scp2(i,j)*1.e-4 ! [m^2]
2831+
garea=scp2(i,j)*A_cgs2mks ! [m^2]
28232832
c
28242833
c --- ----- heat fluxes
28252834
albedo=albs_f*ricclm(i,j,k)+albw(i,j)*(1.-ricclm(i,j,k))
@@ -2838,7 +2847,7 @@ subroutine inifrc_ben02clim
28382847
call xcsum(lht_sum,util3,ip)
28392848
call xcsum(sht_sum,util4,ip)
28402849
c
2841-
fac=1.e4/(12.*area)
2850+
fac=(L_mks2cgs*L_mks2cgs)/(12.*area)
28422851
swa_ave=swa_sum*fac
28432852
lwa_ave=lwa_sum*fac
28442853
lht_ave=lht_sum*fac

‎ben02/sfcstr_ben02.F

+4-3
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
! ------------------------------------------------------------------------------
2-
! Copyright (C) 2004-2020 Mats Bentsen
2+
! Copyright (C) 2004-2022 Mats Bentsen, Mehmet Ilicak
33
!
44
! This file is part of BLOM.
55
!
@@ -24,6 +24,7 @@ subroutine sfcstr_ben02(m,n,mm,nn,k1m,k1n)
2424
c --- ------------------------------------------------------------------
2525
c
2626
use mod_xc
27+
use mod_constants, only: P_mks2cgs
2728
use mod_forcing, only: ztx, mty, taux, tauy
2829
use mod_seaice, only: ficem, hicem, tauxice, tauyice
2930
use mod_checksum, only: csdiag, chksummsk
@@ -44,14 +45,14 @@ subroutine sfcstr_ben02(m,n,mm,nn,k1m,k1n)
4445
do i=max(1,ifu(j,l)),min(ii,ilu(j,l))
4546
facice=(ficem(i,j)+ficem(i-1,j))
4647
. *min(2.,hicem(i,j)+hicem(i-1,j))*.25
47-
taux(i,j)=10.*(ztx(i,j)*(1.-facice)+tauxice(i,j)*facice)
48+
taux(i,j)=P_mks2cgs*(ztx(i,j)*(1.-facice)+tauxice(i,j)*facice)
4849
enddo
4950
enddo
5051
do l=1,isv(j)
5152
do i=max(1,ifv(j,l)),min(ii,ilv(j,l))
5253
facice=(ficem(i,j)+ficem(i,j-1))
5354
. *min(2.,hicem(i,j)+hicem(i,j-1))*.25
54-
tauy(i,j)=10.*(mty(i,j)*(1.-facice)+tauyice(i,j)*facice)
55+
tauy(i,j)=P_mks2cgs*(mty(i,j)*(1.-facice)+tauyice(i,j)*facice)
5556
enddo
5657
enddo
5758
enddo

‎ben02/thermf_ben02.F

+31-26
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
! ------------------------------------------------------------------------------
2-
! Copyright (C) 2002-2021 Mats Bentsen
2+
! Copyright (C) 2002-2022 Mats Bentsen, Mehmet Ilicak
33
!
44
! This file is part of BLOM.
55
!
@@ -21,7 +21,8 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
2121
c
2222
c --- NERSC version of thermf.
2323
c
24-
use mod_constants, only: spcifh, t0deg, epsil, onem
24+
use mod_constants, only: spcifh, t0deg, alpha0, epsilt, onem,
25+
. g2kg, kg2g, L_mks2cgs, M_mks2cgs
2526
use mod_time, only: nday_in_year, nday_of_year, nstep,
2627
. nstep_in_day, baclin,
2728
. xmi, l1mi, l2mi, l3mi, l4mi, l5mi
@@ -71,7 +72,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
7172
. fice,hice,hsnw,tsrf,fice0,hice0,hsnw0,qsww,qnsw,tice,albi,
7273
. tsmlt,albi_h,qswi,dh,qsnwf,fcond,qdamp,qsmlt,qo2i,qbot,swfac,
7374
. dtml,q,volice,df,dvi,dvs,fwflx,sstc,rice,trxflx,sssc,srxflx,
74-
. totsfl,totwfl,sflxc,totsrp,totsrn
75+
. totsfl,totwfl,sflxc,totsrp,totsrn,A_cgs2mks
7576
#ifdef TRC
7677
integer nt
7778
real, dimension(ntr,1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
@@ -82,10 +83,12 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
8283
c
8384
real intp1d
8485
external intp1d
86+
c
87+
A_cgs2mks=1./(L_mks2cgs**2)
8588
c
8689
c --- Due to conservation, the ratio of ice and snow density must be
8790
c --- equal to the ratio of ice and snow heat of fusion
88-
if (abs(fuss/fusi-rhosnw/rhoice).gt.epsil) then
91+
if (abs(fuss/fusi-rhosnw/rhoice).gt.epsilt) then
8992
if (mnproc.eq.1) then
9093
write (lp,*)
9194
. 'thermf: check consistency between snow/ice densities'
@@ -97,7 +100,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
97100
c
98101
c --- Set various constants
99102
dt=baclin ! Time step
100-
cpsw=spcifh*1.e3 ! Specific heat of seawater
103+
cpsw=spcifh*M_mks2cgs ! Specific heat of seawater
101104
rnf_fac=baclin/real(nrfets*86400) ! Runoff reservoar detrainment rate
102105
sag_fac=exp(-sagets*dt) ! Snow aging rate
103106
c
@@ -326,7 +329,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
326329
c --- ----- Ice volume that has to freeze to balance the heat budget
327330
volice=-(qsww+qnsw-q)*(1.-fice)*dt/fusi
328331
c
329-
if (volice.gt.epsil) then
332+
if (volice.gt.epsilt) then
330333
c
331334
c --- ------- New ice in the lead is formed with a specified thickness.
332335
c --- ------- Estimate the change in ice fraction
@@ -344,7 +347,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
344347
c --- ----- If the lead is warming, let the fraction (1 - fice) go to
345348
c --- ----- warm the lead, and the fraction fice to melt ice laterally
346349
fice=fice-(swfac*qsww+qnsw)*fice*dt
347-
. /max(hice*fusi+hsnw*fuss,epsil)
350+
. /max(hice*fusi+hsnw*fuss,epsilt)
348351
if (fice.lt.0.) then
349352
fice=0.
350353
hice=0.
@@ -398,14 +401,14 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
398401
fwflx=eva(i,j)+lip(i,j)+sop(i,j)+rnf(i,j)+rfi(i,j)+fmltfz(i,j)
399402
c
400403
c --- --- Salt flux [kg m-2 s-1] (positive downwards)
401-
sfl(i,j)=-sice*dvi*rhoice/dt*1.e-3
404+
sfl(i,j)=-sice*dvi*rhoice/dt*g2kg
402405
c
403406
c --- --- Salt flux due to brine rejection of freezing sea
404407
c --- --- ice [kg m-2 m-1] (positive downwards)
405-
brnflx(i,j)=max(0.,-sotl*fmltfz(i,j)*1.e-3+sfl(i,j))
408+
brnflx(i,j)=max(0.,-sotl*fmltfz(i,j)*g2kg+sfl(i,j))
406409
c
407410
c --- --- Virtual salt flux [kg m-2 s-1] (positive downwards)
408-
vrtsfl(i,j)=-sotl*fwflx*1.e-3
411+
vrtsfl(i,j)=-sotl*fwflx*g2kg
409412
c
410413
c --- --- Store area weighted virtual salt flux and fresh water flux
411414
util1(i,j)=vrtsfl(i,j)*scp2(i,j)
@@ -415,11 +418,11 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
415418
hmltfz(i,j)=(dvi*fusi+dvs*fuss)/dt
416419
c
417420
c --- --- Total heat flux in BLOM units [W cm-2] (positive upwards)
418-
surflx(i,j)=-(swa(i,j)+nsf(i,j)+hmltfz(i,j))*1.e-4
421+
surflx(i,j)=-(swa(i,j)+nsf(i,j)+hmltfz(i,j))*A_cgs2mks
419422
c
420423
c --- --- Short-wave heat flux in BLOM units [W cm-2] (positive
421424
c --- --- upwards)
422-
sswflx(i,j)=-qsww*(1.-fice0)*1.e-4
425+
sswflx(i,j)=-qsww*(1.-fice0)*A_cgs2mks
423426
c
424427
#ifdef TRC
425428
c --- ------------------------------------------------------------------
@@ -452,7 +455,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
452455
endif
453456
# endif
454457
# endif
455-
trflx(nt,i,j)=-trc(i,j,k1n,nt)*fwflx*1.e-3
458+
trflx(nt,i,j)=-trc(i,j,k1n,nt)*fwflx*g2kg
456459
ttrsf(nt,i,j)=trflx(nt,i,j)*scp2(i,j)
457460
ttrav(nt,i,j)=trc(i,j,k1n,nt)*scp2(i,j)
458461
enddo
@@ -465,16 +468,16 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
465468
surrlx(i,j)=0.
466469
c
467470
c --- --- If trxday>0 , apply relaxation towards observed sst
468-
if (trxday.gt.epsil) then
471+
if (trxday.gt.epsilt) then
469472
sstc=intp1d(sstclm(i,j,l1mi),sstclm(i,j,l2mi),
470473
. sstclm(i,j,l3mi),sstclm(i,j,l4mi),
471474
. sstclm(i,j,l5mi),xmi)
472475
rice=intp1d(ricclm(i,j,l1mi),ricclm(i,j,l2mi),
473476
. ricclm(i,j,l3mi),ricclm(i,j,l4mi),
474477
. ricclm(i,j,l5mi),xmi)
475478
sstc=(1.-rice)*max(sstc,tice_f)+rice*tice_f
476-
trxflx=spcifh*100.*min(hmxl,trxdpt)/(trxday*86400.)
477-
. *min(trxlim,max(-trxlim,sstc-tmxl))
479+
trxflx=spcifh*L_mks2cgs*min(hmxl,trxdpt)/(trxday*86400.)
480+
. *min(trxlim,max(-trxlim,sstc-tmxl))/alpha0
478481
surrlx(i,j)=-trxflx
479482
else
480483
trxflx=0.
@@ -496,12 +499,12 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
496499
salrlx(i,j)=0.
497500
c
498501
c --- --- if srxday>0 , apply relaxation towards observed sss
499-
if (srxday.gt.epsil) then
502+
if (srxday.gt.epsilt) then
500503
sssc=intp1d(sssclm(i,j,l1mi),sssclm(i,j,l2mi),
501504
. sssclm(i,j,l3mi),sssclm(i,j,l4mi),
502505
. sssclm(i,j,l5mi),xmi)
503-
srxflx=100.*min(hmxl,srxdpt)/(srxday*86400.)
504-
. *min(srxlim,max(-srxlim,sssc-smxl))
506+
srxflx=L_mks2cgs*min(hmxl,srxdpt)/(srxday*86400.)
507+
. *min(srxlim,max(-srxlim,sssc-smxl))/alpha0
505508
salrlx(i,j)=-srxflx
506509
util3(i,j)=max(0.,salrlx(i,j))*scp2(i,j)
507510
util4(i,j)=min(0.,salrlx(i,j))*scp2(i,j)
@@ -538,7 +541,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
538541
c --- -------------------------------------------------------------------
539542
c
540543
ustar(i,j)=(min(ustari(i,j),.8e-2)*fice0
541-
. +ustarw(i,j)*(1.-fice0))*1.e2
544+
. +ustarw(i,j)*(1.-fice0))*L_mks2cgs
542545
c
543546
enddo
544547
enddo
@@ -556,7 +559,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
556559
call xcsum(totwfl,util2,ips)
557560
c
558561
c --- Correction for the virtual salt flux [kg m-2 s-1]
559-
sflxc=(-sref*totwfl*1.e-3-totsfl)/area
562+
sflxc=(-sref*totwfl*g2kg-totsfl)/area
560563
if (mnproc.eq.1) then
561564
write (lp,*) 'thermf: totsfl/area,sflxc',totsfl/area,sflxc
562565
endif
@@ -567,8 +570,10 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
567570
do j=1,jj
568571
do l=1,isp(j)
569572
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
570-
salflx(i,j)=-(vrtsfl(i,j)+sflxc+sfl(i,j))*1.e2
571-
brnflx(i,j)=-brnflx(i,j)*1.e2
573+
salflx(i,j)=-(vrtsfl(i,j)+sflxc+sfl(i,j))
574+
. *(kg2g*(M_mks2cgs/L_mks2cgs**2))
575+
brnflx(i,j)=-brnflx(i,j)
576+
. *(kg2g*(M_mks2cgs/L_mks2cgs**2))
572577
enddo
573578
enddo
574579
enddo
@@ -577,7 +582,7 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
577582
c --- if srxday>0 and srxbal=.true. , balance the sss relaxation flux
578583
c --- so the net input of salt in grid cells connected to the world
579584
c --- ocean is zero
580-
if (srxday.gt.epsil.and.srxbal) then
585+
if (srxday.gt.epsilt.and.srxbal) then
581586
call xcsum(totsrp,util3,ipwocn)
582587
call xcsum(totsrn,util4,ipwocn)
583588
if (abs(totsrp).gt.abs(totsrn)) then
@@ -632,14 +637,14 @@ subroutine thermf_ben02(m,n,mm,nn,k1m,k1n)
632637
c tottrav=tottrav/area
633638
cc
634639
c trflxc=(-tottrsf)/area
635-
c trflxc=(-tottrav*totwfl*1.e-3-tottrsf)/area
640+
c trflxc=(-tottrav*totwfl*g2kg-tottrsf)/area
636641
trflxc=0.
637642
c
638643
c$OMP PARALLEL DO PRIVATE(l,i)
639644
do j=1,jj
640645
do l=1,isp(j)
641646
do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
642-
trflx(nt,i,j)=-(trflx(nt,i,j)+trflxc)*1.e2
647+
trflx(nt,i,j)=-(trflx(nt,i,j)+trflxc)*L_mks2cgs
643648
enddo
644649
enddo
645650
enddo

‎cesm/sfcstr_cesm.F

+4-3
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
! ------------------------------------------------------------------------------
2-
! Copyright (C) 2015-2020 Mats Bentsen
2+
! Copyright (C) 2015-2022 Mats Bentsen, Mehmet Ilicak
33
!
44
! This file is part of BLOM.
55
!
@@ -24,6 +24,7 @@ subroutine sfcstr_cesm(m,n,mm,nn,k1m,k1n)
2424
c --- ------------------------------------------------------------------
2525
c
2626
use mod_xc
27+
use mod_constants, only: P_mks2cgs
2728
use mod_forcing, only: ztx, mty, taux, tauy
2829
use mod_checksum, only: csdiag, chksummsk
2930
c
@@ -37,12 +38,12 @@ subroutine sfcstr_cesm(m,n,mm,nn,k1m,k1n)
3738
do j=1,jj
3839
do l=1,isu(j)
3940
do i=max(1,ifu(j,l)),min(ii,ilu(j,l))
40-
taux(i,j)=10.*ztx(i,j)
41+
taux(i,j)=P_mks2cgs*ztx(i,j)
4142
enddo
4243
enddo
4344
do l=1,isv(j)
4445
do i=max(1,ifv(j,l)),min(ii,ilv(j,l))
45-
tauy(i,j)=10.*mty(i,j)
46+
tauy(i,j)=P_mks2cgs*mty(i,j)
4647
enddo
4748
enddo
4849
enddo

0 commit comments

Comments
 (0)