Skip to content

Adding mre terms #242

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
logan-nc merged 39 commits into from
Mar 10, 2026
Merged

Adding mre terms #242

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
39 commits
Select commit Hold shift + click to select a range
3bd70cd
standardising print statements of psihigh, psilow for rdcon and stride
StuartBenjamin Feb 11, 2026
735f304
adding printing of H
StuartBenjamin Feb 11, 2026
fb3e462
adding printing of volume, Zeff, bwall to netcdf
StuartBenjamin Feb 11, 2026
d75fb21
updating printing MRE terms to netcdf
StuartBenjamin Feb 11, 2026
9e5b40a
removing commented out printing of J_bs to netcdf (error in formula)
StuartBenjamin Feb 11, 2026
ed29f9e
adding new variables for MRE terms
StuartBenjamin Feb 11, 2026
05650ce
adding new variables for calculation of local elongation, triangularity
StuartBenjamin Feb 11, 2026
9e6b8b0
adding varying-sized spline for optional inclusion of mre terms
StuartBenjamin Feb 11, 2026
6193d11
initialising extremum surface locations, adding some comments
StuartBenjamin Feb 11, 2026
442c9fe
updating extremum surface locations, adding some comments too
StuartBenjamin Feb 11, 2026
53d1fe0
adding more comments
StuartBenjamin Feb 11, 2026
461253e
adding MRE term computations
StuartBenjamin Feb 11, 2026
120ea66
adding calculation and saving (in mreterms spline) of mre terms, if M…
StuartBenjamin Feb 11, 2026
1e46f43
adding warnings about avg_Jboot_dot_B since formula is incorrect, rem…
StuartBenjamin Feb 11, 2026
4f9ca23
adding qlow, qhigh to rdcon.in for consistency with stride.in
StuartBenjamin Feb 11, 2026
782dd6b
adding new rdcon read in variables Zeff, MRE_flag, geom_flag
StuartBenjamin Feb 11, 2026
db18f63
removing commented out code
StuartBenjamin Feb 11, 2026
c55977f
standardising STRIDE makefile compared to rdcon, to call in $(NETCDF…
StuartBenjamin Feb 11, 2026
f0090e1
adding new integer variables in rdcon_netcdf for MRE terms
StuartBenjamin Feb 11, 2026
c98d017
simplifying/clarifying logic for printing ideal stability information…
StuartBenjamin Feb 11, 2026
1265199
pulling in changes from current develop branch
StuartBenjamin Feb 11, 2026
9f2a524
making variable name (Hbs_prefac) consistent with value being computed
StuartBenjamin Feb 11, 2026
325f04b
making explanation more understandable
StuartBenjamin Feb 11, 2026
47385ea
removing commented out code for Jtor (was incorrect)
StuartBenjamin Feb 11, 2026
7689bb3
avoiding mixing up branches
StuartBenjamin Feb 11, 2026
eb31534
Apply suggestions from code review
StuartBenjamin Feb 11, 2026
341bbe0
Merge remote-tracking branch 'origin/develop' with updated Make scrip…
StuartBenjamin Feb 24, 2026
793c74f
fixing co-pilot's format error
StuartBenjamin Mar 1, 2026
aa5f53b
Adding <|nabla psi|>
StuartBenjamin Mar 1, 2026
a5549d0
removing avg_Jboot_dot_B due to broken formula
StuartBenjamin Mar 1, 2026
91cfd0e
adding nabla_psi integer 'npsi_id'
StuartBenjamin Mar 1, 2026
177924e
printing <|nabla_psi>| to netcdf if geom_flag is true
StuartBenjamin Mar 1, 2026
c929bad
compilation error fix
StuartBenjamin Mar 1, 2026
72b0a0c
Remove automatic Claude PR review workflow
logan-nc Mar 10, 2026
0ce6b6f
Comment out Jtor dead code in mercier_scan
logan-nc Mar 10, 2026
5e92fa7
Fix kind suffix and initial value for extrema arrays in mercier_scan
logan-nc Mar 10, 2026
68120fd
Gate extrema init and tracking on MRE_flag in mercier_scan
logan-nc Mar 10, 2026
d51e7c6
Guard against division by zero at magnetic axis in mercier_scan
logan-nc Mar 10, 2026
571f2b2
Remove unused jbs_id declaration from rdcon_netcdf_out
logan-nc Mar 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 0 additions & 58 deletions .github/workflows/claude-code-review.yml
View file
Open in desktop

This file was deleted.

124 changes: 91 additions & 33 deletions rdcon/mercier.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,17 +29,20 @@ MODULE rdcon_mercier_mod
SUBROUTINE mercier_scan

INTEGER :: ipsi,itheta
REAL(r8) :: bsq,chi1,di,dpsisq,eta,h,jac,p1,psifac,q,q1,r,
REAL(r8) :: bsq,chi1,di,dpsisq,eta,h,jac,p1,psifac,q,q1,r,z,
$ rfac,term,theta,twopif,v1,v2,v21,v22,v23,v33,
$ bt,f1,Dnc,Dnc_prefac,Wc_prefac,eps_loc,Jpara,ftr,M,
$ Hbs,avg_Jboot_dot_B,mufrac,taua_prefac,taur_prefac
$ Hbs_prefac,mufrac,taua_prefac,taur_prefac,
$ rmean_loc,amean_loc,deltatop_loc,deltabot_loc,delta_loc,
$ eps_eff
REAL(r8), DIMENSION(:), POINTER :: avg
REAL(r8), DIMENSION(2) :: rmax_loc, rmin_loc, zmax_loc, zmin_loc
TYPE(spline_type), TARGET :: ff
c-----------------------------------------------------------------------
c prepare spline types.
c-----------------------------------------------------------------------
IF(MRE_flag)THEN
CALL spline_alloc(ff,mtheta,21)
CALL spline_alloc(ff,mtheta,22)
ELSE
CALL spline_alloc(ff,mtheta,5)
ENDIF
Expand All @@ -58,6 +61,15 @@ SUBROUTINE mercier_scan
q1=sq%fs1(ipsi,4)
chi1=twopi*psio ! d(poloidal flux)/dpsi, chi=poloidal flux
c-----------------------------------------------------------------------
c initialise extremum surface locations (MRE terms only).
c-----------------------------------------------------------------------
IF(MRE_flag)THEN
rmax_loc=-HUGE(1.0_r8)
rmin_loc=HUGE(1.0_r8)
zmax_loc=-HUGE(1.0_r8)
zmin_loc=HUGE(1.0_r8)
ENDIF
c-----------------------------------------------------------------------
c evaluate coordinates and jacobian.
c-----------------------------------------------------------------------
DO itheta=0,mtheta
Expand All @@ -66,9 +78,31 @@ SUBROUTINE mercier_scan
rfac=SQRT(rzphi%f(1)) ! minor radius
eta=twopi*(theta+rzphi%f(2)) ! machine poloidal angle
r=ro+rfac*COS(eta) ! major radius R
z=zo+rfac*SIN(eta) ! vertical coordinate Z
jac=rzphi%f(4) ! jacobian of mag. coordinates
bt=twopif/(twopi*r) ! toroidal B field
c-----------------------------------------------------------------------
c update extremum surface locations (MRE terms only).
c-----------------------------------------------------------------------
IF(MRE_flag)THEN
IF(r>rmax_loc(1))THEN
rmax_loc(1)=r
rmax_loc(2)=z
ENDIF
IF(r<rmin_loc(1))THEN
rmin_loc(1)=r
rmin_loc(2)=z
ENDIF
IF(z>zmax_loc(2))THEN
zmax_loc(1)=r
zmax_loc(2)=z
ENDIF
IF(z<zmin_loc(2))THEN
zmin_loc(1)=r
zmin_loc(2)=z
ENDIF
ENDIF
c-----------------------------------------------------------------------
c evaluate other local quantities.
c-----------------------------------------------------------------------
v21=rzphi%fy(1)/(2*rfac*jac)
Expand Down Expand Up @@ -102,6 +136,7 @@ SUBROUTINE mercier_scan
ff%fs(itheta,19)=SQRT(dpsisq)/r
ff%fs(itheta,20)=r**2*v1/jac !overbar{R^2} (Hegna 1999)
ff%fs(itheta,21)=r**2 !avg{R^2} ~ [m^2] (Hegna 1999)
ff%fs(itheta,22)=SQRT(dpsisq)
ENDIF
ff%fs(itheta,:)=ff%fs(itheta,:)*jac/v1
ENDDO
Expand Down Expand Up @@ -139,21 +174,41 @@ SUBROUTINE mercier_scan
taur_prefac=avg(1)/avg(5)*mu0
!to get taur, divide by local resistivity (see resist.f)
c-----------------------------------------------------------------------
c compute surface geometric quantities
c-----------------------------------------------------------------------
rmean_loc=0.5d0*(rmax_loc(1)+rmin_loc(1))
amean_loc=0.5d0*(rmax_loc(1)-rmin_loc(1))
eps_loc=amean_loc/rmean_loc
c at the magnetic axis amean_loc=0: triangularity is undefined
c (->inf) and eps_loc=0 so eps_eff=0 and ftr=0 regardless of delta
IF(amean_loc > 0.0_r8)THEN
deltatop_loc=(rmean_loc-zmax_loc(1))/amean_loc
deltabot_loc=(rmean_loc-zmin_loc(1))/amean_loc
delta_loc=0.5d0*(deltatop_loc+deltabot_loc)
eps_eff=MAX(0.d0,
$ 0.67d0*(1.d0-1.4d0*delta_loc*ABS(delta_loc))*eps_loc)
ELSE
deltatop_loc=HUGE(1.0_r8)
deltabot_loc=HUGE(1.0_r8)
delta_loc=HUGE(1.0_r8)
eps_eff=0.0_r8
ENDIF
c-----------------------------------------------------------------------
c simple estimates of trapped fraction from Sauter 2002:
c <https://infoscience.epfl.ch/server/api/core/bitstreams/c42baba0-9
c 909-4f21-978a-0d0c2646c3ad/content>
c-----------------------------------------------------------------------
eps_loc=avg(11)/avg(12)
ftr=1.d0-(1-eps_loc)**2/
$ (SQRT(1-eps_loc**2)*(1.d0+1.46d0*SQRT(eps_loc)))
ftr=1.d0-SQRT((1.d0-eps_loc)/(1.d0+eps_loc))*(1.d0-eps_eff)/
$ (1.d0+2.d0*SQRT(eps_eff))
ftr=MIN(ftr,1.0d0)
c-----------------------------------------------------------------------
c simple estimates of Jboot and bootstrap drive from
c Callen, 2010 UW-CPTC 09-6R, and Hegna 1999
c simple estimates of bootstrap drive-related terms from
c Callen, 2010 UW-CPTC 09-6R, and Hegna PoP 1999
c-----------------------------------------------------------------------
! mu_e_on_nu_e=(ftr/(1.d0-ftr))*(1.d0+0.533d0/Zeff)
! ^ Taking banana limit of eq. B17 (& B14) in Callen
mufrac=ftr*(1.d0+0.533d0/Zeff)/
$ ((1.d0-ftr)+ftr*(1.d0+0.533d0/Zeff))
avg_Jboot_dot_B=-mufrac*(twopif/chi1)*p1 !mu0 included in p1
c-----------------------------------------------------------------------
c evaluate geometric prefactors of MRE stability terms from
c Hegna 1999 https://doi.org/10.1063/1.873661
Expand All @@ -166,32 +221,35 @@ SUBROUTINE mercier_scan
c evaluate geometric prefactor of Wc from
c Schlutt and Hegna 2012 https://doi.org/10.1063/1.4747500
c-----------------------------------------------------------------------
Wc_prefac=(v1*avg(5)/(q*psio))*
$ v1*avg(7)*(q*psio)*
$ q**6/((q1/psio)**2)
!Divide by mode num. m^2 to finish, units are Wb^4
!Note, psi coordinate in original paper is toroidal flux
Wc_prefac=(v1*avg(5)/(q*psio))*! overbar{J B^2} [in psi_tor]
$ v1*avg(7)*(q*psio)* ! overbar{J |nabla psi|^2} [in psi_tor]
$ (q**6/((q1/psio)**2))* ! {d iota / d psi_tor}^(-2)[in psi_tor]
$ (1/(q*psio))**4 ! *Conversion from psi_tor^4 to psi_pol_norm^4
!Divide by mode num. m^2 to finish, units are (psi_pol_norm)^4
!Note, psi coordinate in original paper is toroidal flux
!per radian `psi_tor'. *Using small island approximation.
c-----------------------------------------------------------------------
c parallel current density from Freidberg Ideal MHD eqs. 6.15, 6.16.
c note psi in eqs. 6.15, 6.16 is poloidal flux/(2pi), same as psi_in
c-----------------------------------------------------------------------
Jpara=psio*f1*avg(16) + p1*avg(18)*twopif/(twopi*psio) +
Jpara=psio*f1*avg(16)+p1*avg(18)*twopif/(twopi*psio) +
$ (twopif/twopi)**2*f1*avg(17)/psio !(mu0 included in p1))
c Jtor=p1*avg(12)/psio + (twopif/twopi)*f1*avg(13)/psio
c-----------------------------------------------------------------------
c compute H_bs defined in Shi et al. 2024, using identity from
c compute Hbs_prefac from Shi et al. 2024, using identity from
c Glasser et al. 1975. Shi doi -> https://doi.org/10.1063/5.0183474
c-----------------------------------------------------------------------
Hbs=(avg(1)/avg(5))* ! [m^2]
$ (-v1/(twopi**2*psio**2*q1))*avg_Jboot_dot_B ! [m^-5*m^3]
Hbs_prefac=(avg(1)/avg(5))*
$ (-v1/(twopi**2*psio**2*q1))! mult by avg_Jboot_dot_B to get Hbs
c-----------------------------------------------------------------------
c save terms.
c-----------------------------------------------------------------------
mreterms%fs(ipsi,1)=Hbs
mreterms%fs(ipsi,1)=Hbs_prefac
mreterms%fs(ipsi,2)=taua_prefac
mreterms%fs(ipsi,3)=taur_prefac
mreterms%fs(ipsi,4)=ftr
mreterms%fs(ipsi,5)=mufrac
mreterms%fs(ipsi,6)=avg_Jboot_dot_B
mreterms%fs(ipsi,6)=avg(22) ! <|nabla psi|>
mreterms%fs(ipsi,7)=Dnc
mreterms%fs(ipsi,8)=Wc_prefac
mreterms%fs(ipsi,9)=Jpara
Expand All @@ -204,19 +262,19 @@ SUBROUTINE mercier_scan
mreterms%fs(ipsi,16)=avg(20) !overbar{R^2} (Hegna 1999)
mreterms%fs(ipsi,17)=avg(21) !avg{R^2} ~ [m^2] (Hegna 1999)
! will only print out the following if geom_flag is true:
mreterms%fs(ipsi,18)=avg(1)
mreterms%fs(ipsi,19)=avg(2)
mreterms%fs(ipsi,20)=avg(3)
mreterms%fs(ipsi,21)=avg(4)
mreterms%fs(ipsi,22)=avg(5)
mreterms%fs(ipsi,23)=avg(6)
mreterms%fs(ipsi,24)=avg(7)
mreterms%fs(ipsi,25)=avg(14)
mreterms%fs(ipsi,26)=avg(15)
mreterms%fs(ipsi,27)=avg(16)
mreterms%fs(ipsi,28)=avg(17)
mreterms%fs(ipsi,29)=avg(18)
mreterms%fs(ipsi,30)=avg(19)
mreterms%fs(ipsi,18)=avg(1) ! <bsq/dpsisq>
mreterms%fs(ipsi,19)=avg(2) ! <1/dpsisq>
mreterms%fs(ipsi,20)=avg(3) ! <1/bsq>
mreterms%fs(ipsi,21)=avg(4) ! <1/(bsq*dpsisq)>
mreterms%fs(ipsi,22)=avg(5) ! <bsq>
mreterms%fs(ipsi,23)=avg(6) ! <dpsisq/bsq>
mreterms%fs(ipsi,24)=avg(7) ! <dpsisq>
mreterms%fs(ipsi,25)=avg(14)! <dpsisq/(R^2)>
mreterms%fs(ipsi,26)=avg(15)! <1/(R^2)>
mreterms%fs(ipsi,27)=avg(16)! <dpsisq/(R^2 B))>
mreterms%fs(ipsi,28)=avg(17)! <1/(R^2 B)>
mreterms%fs(ipsi,29)=avg(18)! <1/B>
mreterms%fs(ipsi,30)=avg(19)! <sqrt(dpsisq)/R>
ENDIF
120 FORMAT(19(E30.15,1X))
ENDDO
Expand Down
22 changes: 13 additions & 9 deletions rdcon/rdcon_netcdf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,16 +57,15 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)

INTEGER :: i, ncid,
$ i_dim, m_dim, mo_dim, p_dim, i_id, m_id, mo_id, p_id,
$ f_id, q_id, dv_id, mu_id, di_id, dr_id, ca_id,
$ f_id, q_id, dv_id, mu_id, di_id, dr_id, h_id, ca_id,
$ wp_id, wpv_id, wv_id, wvv_id, wt_id, wtv_id, wt0_id,
$ l_dim, l_id, coil_dim, coil_id, dpc_id, dc_id,
$ lp_dim, lp_id, r_dim, r_id, rp_dim, rp_id, pr_id, qr_id,
$ dp_id, ap_id, bp_id, gp_id, dpp_id, lrc_dim, lrc_id
INTEGER :: hbs_id, ta_id, tr_id, ftr_id, muf_id, jbs_id, dnc_id,
INTEGER :: hbs_id, ta_id, tr_id, ftr_id, muf_id, dnc_id,
$ wc_id, jp_id, b_id, bt_id, bpo_id, mir_id, mar_id, mair_id,
$ obr_id, ars_id, a1_id, a2_id, a3_id, a4_id, a5_id,
$ a6_id, a7_id, a14_id, a15_id, a16_id, a17_id, a18_id, a19_id

$ obr_id, ars_id, a1_id, a2_id, a3_id, a4_id, a5_id, a6_id,
$ a7_id, a14_id, a15_id, a16_id, a17_id, a18_id, a19_id, npsi_id
REAL(r4) :: cpusec, wallsec
CHARACTER(2) :: sn
CHARACTER(64) :: ncfile
Expand Down Expand Up @@ -129,6 +128,9 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
CALL check( nf90_put_att(ncid,nf90_global,'betap3', betap3))
CALL check( nf90_put_att(ncid,nf90_global,'betat', betat))
CALL check( nf90_put_att(ncid,nf90_global,'betan', betan))
CALL check( nf90_put_att(ncid,nf90_global,'volume', volume))
CALL check( nf90_put_att(ncid,nf90_global,'Zeff', Zeff))
CALL check( nf90_put_att(ncid,nf90_global,'bwall', bwall))
Comment on lines +131 to +133
Copy link

Copilot AI Feb 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Global NetCDF attributes in this file consistently use lowercase names (e.g., betan, bwall, volume), but this adds 'Zeff' with mixed case. Consider renaming the attribute to zeff (or otherwise matching the established casing) to avoid surprising downstream readers that treat attribute names as case-sensitive.

Copilot uses AI. Check for mistakes.
Copy link
Contributor

@logan-nc logan-nc Mar 10, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The capitalization is ok here. It follows standard annotations in the field

CALL check( nf90_put_att(ncid,nf90_global,'bt0', bt0))
CALL check( nf90_put_att(ncid,nf90_global,'q0', q0))
CALL check( nf90_put_att(ncid,nf90_global,'qmin', qmin))
Expand Down Expand Up @@ -217,10 +219,11 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
$ "Safety Factor") )
CALL check( nf90_def_var(ncid, "di", nf90_double, p_dim, di_id) )
CALL check( nf90_def_var(ncid, "dr", nf90_double, p_dim, dr_id) )
CALL check( nf90_def_var(ncid, "h", nf90_double, p_dim, h_id) )
CALL check( nf90_def_var(ncid, "ca1", nf90_double, p_dim, ca_id))
IF(MRE_flag)THEN
CALL check( nf90_def_var(ncid,
$ "Hbs", nf90_double, p_dim, hbs_id))
$ "Hbs_prefac", nf90_double, p_dim, hbs_id))
CALL check( nf90_def_var(ncid,
$ "tau_a", nf90_double, p_dim, ta_id))
CALL check( nf90_def_var(ncid,
Expand All @@ -230,8 +233,6 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
CALL check( nf90_def_var(ncid,
$ "mufrac", nf90_double, p_dim, muf_id))
CALL check( nf90_def_var(ncid,
$ "avg_mu0Jbs_dot_B", nf90_double, p_dim, jbs_id))
CALL check( nf90_def_var(ncid,
$ "Dnc", nf90_double, p_dim, dnc_id))
CALL check( nf90_def_var(ncid,
$ "Wc", nf90_double, p_dim, wc_id))
Copy link

Copilot AI Feb 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

avg_mu0Jbs_dot_B is no longer defined/written, but jbs_id is still declared (and mreterms%fs(:,6) still exists). Consider removing jbs_id (and/or documenting the skipped index) so the netCDF writer and mreterms layout don’t silently diverge.

Suggested change
$ "Wc", nf90_double, p_dim, wc_id))
$ "Wc", nf90_double, p_dim, wc_id))
c
c NOTE: The MRE terms layout (e.g., mreterms%fs(:,6)) still reserves
c an index for the deprecated quantity avg_mu0Jbs_dot_B. That term
c is intentionally no longer written to this NetCDF file, so no
c nf90_def_var call appears here for it. The sequence of variables
c defined below therefore skips that index by design.
c

Copilot uses AI. Check for mistakes.
Copy link
Contributor

@logan-nc logan-nc Mar 10, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is overly verbose. We should just remove the index from being declared and add a one-line annotation wherever mreterms%fs(:,6) is actually calculated saying it comes from [citation] but gives suspect results and is thus not output for now

Copy link
Collaborator Author

@StuartBenjamin StuartBenjamin Mar 10, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This copilot comment is out of date and should be ignored -> mreterms%fs(:,6) now stores avg(|nabla psi|) which is correctly calculated, and printed on lines 284, 285, and 417

Expand Down Expand Up @@ -280,6 +281,8 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
$ "avg_18", nf90_double, p_dim, a18_id))
CALL check( nf90_def_var(ncid,
$ "avg_19", nf90_double, p_dim, a19_id))
CALL check( nf90_def_var(ncid,
$ "avg_nabla_psi", nf90_double, p_dim, npsi_id))
ENDIF
ENDIF
IF(ode_flag .AND. vac_flag)THEN !shift to .OR.
Expand Down Expand Up @@ -377,6 +380,7 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
CALL check( nf90_put_var(ncid,q_id, sq%fs(:,4)))
CALL check( nf90_put_var(ncid,di_id, locstab%fs(:,1)/sq%xs(:)))
CALL check( nf90_put_var(ncid,dr_id, locstab%fs(:,2)/sq%xs(:)))
CALL check( nf90_put_var(ncid,h_id, locstab%fs(:,3)))
CALL check( nf90_put_var(ncid,ca_id, locstab%fs(:,4)))

IF(MRE_flag)THEN
Expand All @@ -385,7 +389,6 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
CALL check( nf90_put_var(ncid,tr_id, mreterms%fs(:,3)))
CALL check( nf90_put_var(ncid,ftr_id, mreterms%fs(:,4)))
CALL check( nf90_put_var(ncid,muf_id, mreterms%fs(:,5)))
CALL check( nf90_put_var(ncid,jbs_id, mreterms%fs(:,6)))
CALL check( nf90_put_var(ncid,dnc_id, mreterms%fs(:,7)))
CALL check( nf90_put_var(ncid,wc_id, mreterms%fs(:,8)))
CALL check( nf90_put_var(ncid,jp_id, mreterms%fs(:,9)))
Expand All @@ -411,6 +414,7 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
CALL check( nf90_put_var(ncid,a17_id, mreterms%fs(:,28)))
CALL check( nf90_put_var(ncid,a18_id, mreterms%fs(:,29)))
CALL check( nf90_put_var(ncid,a19_id, mreterms%fs(:,30)))
CALL check( nf90_put_var(ncid,npsi_id, mreterms%fs(:,6)))
ENDIF
ENDIF

Expand Down
Loading