Opened 3 years ago
Closed 3 years ago
#2669 closed Bug (fixed)
ISF: Incorrect ocean temperature trend when ice-shelf melting in the trunk
| Reported by: | clem | Owned by: | mathiot |
|---|---|---|---|
| Priority: | normal | Milestone: | |
| Component: | ICB | Version: | trunk |
| Severity: | minor | Keywords: | |
| Cc: |
Description
The ocean temperature trend coming from isf melting is wrong when cavities are closed and melting is parameterized. The part that is wrong is the sensible flux, so I guess it was undetected since it is small wrt latent heat flux. Moreover the way this trend is written is extremely confusing. In the following, I try to explain the best I can.
If I understand correctly, the ocean temperature "trend" must be written as: fwf * (cp*T-Lfus) / (rho*cp) in degC*m/s, with fwf being the positive freshwater, T the freezing temperature and Lfus the latent heat of fusion. Hence melting is responsible for cooling the surrounding ocean.
In ISF routines, there are two ways of defining temperature "trend": in isfcav.F90 if the cavity is open or in isfpar.F90 if the cavity is closed and melting is parameterized. Sign is opposite in the two routines and if I did my maths ok, I end up with temperature trends as:
fwf * (cp*T-Lfus) / (rho*cp) for isfcav => ok
fwf * (-cp*T-Lfus) / (rho*cp) for isfpar => not ok (sensible flux has the wrong sign)
In details, in isfcav.F90 the temperature trend is defined as:
zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) ) ... ! set temperature content ptsc(:,:,jp_tem) = -zqh(:,:) * r1_rho0_rcp
with qhc (sensible flux) and qoce (latent flux) defined in isfcavmlt.F90 (here for specified melt):
pqfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1) ! fwf ( >0 out) pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocean heat flux ( >0 out) pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( >0 out)
In isfpar.F90 the temperature trend is defined as:
zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) ) ... ! set temperature content ptsc(:,:,jp_tem) = zqh(:,:) * r1_rho0_rcp
with qhc and qoce defined in isfparmlt.F90 (here for specified melt):
pqfwf(:,:) = - sf_isfpar_fwf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) pqoce(:,:) = pqfwf(:,:) * rLfusisf ! ocean/ice shelf flux assume to be equal to latent heat flux pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux
I see a couple of other problems:
1) Signs are confusing (not the same between isfpar and isfcav).
2) The freshwater flux (qfwf or fwfisf) is negative when melting. It does not really follow the nemo general rule of having positive runoffs.
Fix
Correcting temperature trend and solving point 1) is pretty easy. Solving point 2) is much harder since isf freshwater is used in many places of the code and we do not have configurations to test all of the options.
I recommend to first correct the temperature trend, then solve the signs confusion and maybe change the sign of freshwater flux
Commit History (11)
| Changeset | Author | Time | ChangeLog |
|---|---|---|---|
| 14995 | mathiot | 2021-06-15T19:15:26+02:00 | ticket #2669 : merge ticket2669_isf_fluxes into trunk |
| 14994 | mathiot | 2021-06-15T16:39:31+02:00 | ticket #2669: update to the head of trunk |
| 14989 | mathiot | 2021-06-14T15:30:07+02:00 | ticket #2669: add extra comments in the headers |
| 14988 | mathiot | 2021-06-14T15:20:22+02:00 | ticket #2669: isf outputs are still in the SBC part of field_def and it should not as ISF no more in surface module + change long_name of isf flux (less confusing) |
| 14926 | mathiot | 2021-05-28T17:22:05+02:00 | ticket #2669 : change long_name att of fwfisf* variables |
| 14924 | mathiot | 2021-05-28T16:07:37+02:00 | ticket #2669 : off topic again, I plugged in the fix for the oasis case made by Christrian. |
| 14923 | mathiot | 2021-05-28T15:43:29+02:00 | ticket #2669: a bit off topic, but I updated the ISOMIP+ MY_SRC to match the trunk (ease comparison of sources) |
| 14919 | mathiot | 2021-05-27T20:01:49+02:00 | ticket #2669: set sign convention of thermald output in bg03 the same as in other mlt scheme |
| 14917 | mathiot | 2021-05-27T19:30:09+02:00 | ticket #2669: add masking of isftfrz_* variables, fix issue on vertical mean in bg03 case, fix sign of qlatisf output |
| 14916 | mathiot | 2021-05-27T18:40:14+02:00 | ticket #2669: bg03 option uses rn_gammat0 defined in cav namelist bloc (corrected by adding a new namelist parameter rn_isfpar_bg03_gt0) |
| 14904 | mathiot | 2021-05-26T17:26:17+02:00 | creation of branch to work on ticket #2669 |
Change History (22)
comment:1 Changed 3 years ago by clem
comment:2 Changed 3 years ago by clem
Well, for bg03 one must read:
actual formulation:
pqoce(:,:) = rho0 * rcp * rn_gammat0 * risfLeff(:,:) * e1t(:,:) * ( ztavg(:,:) - ztfrz(:,:) ) * r1_e1e2t(:,:) pqfwf(:,:) = - pqoce(:,:) / rLfusisf ! derived from the latent heat flux pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux
instead of
pqoce(:,:) = - rho0 * rcp * rn_gammat0 * risfLeff(:,:) * e1t(:,:) * ( ztavg(:,:) - ztfrz(:,:) ) * r1_e1e2t(:,:) pqfwf(:,:) = pqoce(:,:) / rLfusisf ! derived from the latent heat flux pqhc (:,:) = - pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux
comment:3 Changed 3 years ago by mathiot
- Owner set to mathiot
- Status changed from new to assigned
comment:4 Changed 3 years ago by mathiot
I will look at the 'oasis' one but no test will be made to test it for open and closed cavities. As mentioned in the doc:
1210 The \forcode{'oasis'} is a prototype of what could be a method to spread precipitation on Antarctic ice sheet as ice shelf melt inside the cavity when a coupled model Atmosphere/Ocean? is used.
1211 It has not been tested and therefore the model will stop if you try to use it.
For the 'bg03', this bits of code is from DRAKKAR project in ~2009 and included in the ice shelf module, but has not been re-tested since ~2014 (I think). Based on recent tests by Favier et al. (2019) and work from C. Burgard, with N. Jourdain we do not advise it at all (a quadratic param give much better results).
In the doc it is mentioned like this:
1338 This parametrisation has not been tested since a while and based on \citet{Favier2019},
1339 this parametrisation should probably not be used.
I will look at it but I will not run any test. To be be more clear, I will add a stop in the stop to prevent its usage and mentioned it is a parametrisation exemple which could be useful as starting point.
I will check the other options with WED025.
comment:5 Changed 3 years ago by mathiot
- I agree with 'the way this trend is written is extremely confusing' from the ticket description, I will correct it. Moreover not doing it, implies a different convention in flux output for 'par' and 'cav'
- I agree with the sign issue of the sensible heat for 'par' 'spe' case. I will run a test with uniform temperature (T0), latent heat switch off and zfrz set to T0 to see if T change or not. It should not if sign is correct.
- I agree with the sign issue of bg03. In the actual formulation melt is > 0 and it should be < 0.
- bg03 flux definition will slightly be changed to have the same order of flux definition as other scheme (fwf, then qoce then qhc). It eases reading and check of sign as at the end, signe for qoce and qhc have to be the same for all the scheme.
Extra changes:
- I will also add some possible debug print for 'par', in order to have the same for 'cav' and 'par'
-
src/OCE/ISF/isfparmlt.F90
9 9 10 10 USE isf_oce ! ice shelf 11 11 USE isftbl , ONLY: isf_tbl ! ice shelf depth average 12 USE isfutils,ONLY: debug ! debug subroutine 12 13 13 14 USE dom_oce ! ocean space and time domain 14 15 USE oce , ONLY: ts ! ocean dynamics and tracers … … 38 39 ! -------------------------------- PUBLIC SUBROUTINE ---------------------------------------------------- 39 40 ! ------------------------------------------------------------------------------------------------------- 40 41 41 SUBROUTINE isfpar_mlt( kt, Kmm, pqhc, pqoce, pqfwf )42 SUBROUTINE isfpar_mlt( kt, Kmm, pqhc, pqoce, pqfwf ) 42 43 !!--------------------------------------------------------------------- 43 44 !! *** ROUTINE isfpar_mlt *** 44 45 !! … … 68 69 CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfpar (should not see this)') 69 70 END SELECT 70 71 ! 72 IF (ln_isfdebug) THEN 73 IF(lwp) WRITE(numout,*) '' 74 CALL debug( 'isfpar_mlt qhc :', pqhc (:,:) ) 75 CALL debug( 'isfpar_mlt qoce :', pqoce(:,:) ) 76 CALL debug( 'isfpar_mlt qfwf :', pqfwf(:,:) ) 77 IF(lwp) WRITE(numout,*) '' 78 END IF 79 ! 71 80 END SUBROUTINE isfpar_mlt
- sign convention will be added to field_def_oce of isf fluxes in long_name attribute.
-
cfgs/SHARED/field_def_nemo-oce.xml
399 399 <field id="ssh_ib" long_name="Inverse barometer sea surface height" standard_name="sea_surface_height_correction_due_to_air_pressure_at_low_frequency" unit="m" /> 400 400 401 401 <!-- * variable related to ice shelf forcing * --> 402 403 <!-- * fwf * --> 404 <field id="fwfisf_cav" long_name="Ice shelf melt rate ( > 0 out )" unit="kg/m2/s" /> 405 <field id="fwfisf_par" long_name="Ice shelf melt rate ( > 0 out )" unit="kg/m2/s" /> 406 <field id="fwfisf3d_cav" long_name="Ice shelf melt rate ( > 0 out )" unit="kg/m2/s" grid_ref="grid_T_3D" /> 407 <field id="fwfisf3d_par" long_name="Ice shelf melt rate ( > 0 out )" unit="kg/m2/s" grid_ref="grid_T_3D" /> 408 409 <!-- * heat fluxes * --> 410 <field id="qoceisf_cav" long_name="Ice shelf ocean heat flux ( > 0 out )" unit="W/m2" /> 411 <field id="qoceisf_par" long_name="Ice shelf ocean heat flux ( > 0 out )" unit="W/m2" /> 412 <field id="qlatisf_cav" long_name="Ice shelf latent heat flux ( > 0 out )" unit="W/m2" /> 413 <field id="qlatisf_par" long_name="Ice shelf latent heat flux ( > 0 out )" unit="W/m2" /> 414 <field id="qhcisf_cav" long_name="Ice shelf heat content flux of injected water ( > 0 out )" unit="W/m2" /> 415 <field id="qhcisf_par" long_name="Ice shelf heat content flux of injected water ( > 0 out )" unit="W/m2" /> 416 <field id="qoceisf3d_cav" long_name="Ice shelf ocean heat flux ( > 0 out )" unit="W/m2" grid_ref="grid_T_3D" /> 417 <field id="qoceisf3d_par" long_name="Ice shelf ocean heat flux ( > 0 out )" unit="W/m2" grid_ref="grid_T_3D" /> 418 <field id="qlatisf3d_cav" long_name="Ice shelf latent heat flux ( > 0 out )" unit="W/m2" grid_ref="grid_T_3D" /> 419 <field id="qlatisf3d_par" long_name="Ice shelf latent heat flux ( > 0 out )" unit="W/m2" grid_ref="grid_T_3D" /> 420 <field id="qhcisf3d_cav" long_name="Ice shelf heat content flux of injected water ( > 0 out )" unit="W/m2" grid_ref="grid_T_3D" /> 421 <field id="qhcisf3d_par" long_name="Ice shelf heat content flux of injected water ( > 0 out )" unit="W/m2" grid_ref="grid_T_3D" /> 422 <field id="qconisf" long_name="Conductive heat flux through the ice shelf ( > 0 out )" unit="W/m2" /> 423 424 <!-- top boundary layer properties --> 402 425 <field id="isftfrz_cav" long_name="freezing point temperature at ocean/isf interface" unit="degC" /> 403 426 <field id="isftfrz_par" long_name="freezing point temperature in the parametrization boundary layer" unit="degC" /> 404 <field id="fwfisf_cav" long_name="Ice shelf melt rate" unit="kg/m2/s" />405 <field id="fwfisf_par" long_name="Ice shelf melt rate" unit="kg/m2/s" />406 <field id="qoceisf_cav" long_name="Ice shelf ocean heat flux" unit="W/m2" />407 <field id="qoceisf_par" long_name="Ice shelf ocean heat flux" unit="W/m2" />408 <field id="qlatisf_cav" long_name="Ice shelf latent heat flux" unit="W/m2" />409 <field id="qlatisf_par" long_name="Ice shelf latent heat flux" unit="W/m2" />410 <field id="qhcisf_cav" long_name="Ice shelf heat content flux of injected water" unit="W/m2" />411 <field id="qhcisf_par" long_name="Ice shelf heat content flux of injected water" unit="W/m2" />412 <field id="fwfisf3d_cav" long_name="Ice shelf melt rate" unit="kg/m2/s" grid_ref="grid_T_3D" />413 <field id="fwfisf3d_par" long_name="Ice shelf melt rate" unit="kg/m2/s" grid_ref="grid_T_3D" />414 <field id="qoceisf3d_cav" long_name="Ice shelf ocean heat flux" unit="W/m2" grid_ref="grid_T_3D" />415 <field id="qoceisf3d_par" long_name="Ice shelf ocean heat flux" unit="W/m2" grid_ref="grid_T_3D" />416 <field id="qlatisf3d_cav" long_name="Ice shelf latent heat flux" unit="W/m2" grid_ref="grid_T_3D" />417 <field id="qlatisf3d_par" long_name="Ice shelf latent heat flux" unit="W/m2" grid_ref="grid_T_3D" />418 <field id="qhcisf3d_cav" long_name="Ice shelf heat content flux of injected water" unit="W/m2" grid_ref="grid_T_3D" />419 <field id="qhcisf3d_par" long_name="Ice shelf heat content flux of injected water" unit="W/m2" grid_ref="grid_T_3D" />420 <field id="ttbl_cav" long_name="temperature in Losch tbl" unit="degC" />421 <field id="ttbl_par" long_name="temperature in the parametrisation boundary layer" unit="degC" />422 427 <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting" unit="degC" /> 423 428 <field id="isfthermald_par" long_name="thermal driving of ice shelf melting" unit="degC" /> 424 429 <field id="isfgammat" long_name="Ice shelf heat-transfert velocity" unit="m/s" /> 425 430 <field id="isfgammas" long_name="Ice shelf salt-transfert velocity" unit="m/s" /> 431 <field id="ttbl_cav" long_name="temperature in Losch tbl" unit="degC" /> 432 <field id="ttbl_par" long_name="temperature in the parametrisation boundary layer" unit="degC" /> 426 433 <field id="stbl" long_name="salinity in the Losh tbl" unit="1e-3" /> 427 434 <field id="utbl" long_name="zonal current in the Losh tbl at T point" unit="m/s" /> 428 435 <field id="vtbl" long_name="merid current in the Losh tbl at T point" unit="m/s" /> 429 436 <field id="isfustar" long_name="ustar at T point used in ice shelf melting" unit="m/s" /> 430 <field id="qconisf" long_name="Conductive heat flux through the ice shelf" unit="W/m2" />431 437 432 438 <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> 433 439 <field id="rho_air" long_name="Air density at 10m above sea surface" standard_name="rho_air_10m" unit="kg/m3" />
comment:6 Changed 3 years ago by mathiot
In 14904:
comment:7 Changed 3 years ago by mathiot
comment:8 Changed 3 years ago by mathiot
I forgot to add the ticket number in the messages:
- commit revision r14906 : update long_name of isf flux variable (add sign convention) and sign convention of the isf melt forcing in namelist
- commit revision r14907 : add some debug diagnostics and improve readability of debug print
- commit revision r14908 : fix consistency between cav and par cases and sign issues
comment:9 Changed 3 years ago by mathiot
In 14916:
comment:10 Changed 3 years ago by mathiot
In 14917:
comment:11 Changed 3 years ago by mathiot
In 14919:
comment:12 Changed 3 years ago by mathiot
In 14923:
comment:13 Changed 3 years ago by mathiot
In 14924:
comment:14 Changed 3 years ago by mathiot
Changes made for this ticket passes sette test and only results for WED025 changed compare to trunk as expected (heat content flux sign):
Current code is : NEMO/branches/2021/ticket2669_isf_fluxes @ r14904 ( last change @ r14904 )
SETTE validation report generated for :
NEMO/branches/2021/ticket2669_isf_fluxes @ r14904 (last changed revision)
on occigen arch file
!!---------------1st pass------------------!!
!----restart----!
WGYRE_PISCES_ST run.stat restartability passed : 14904
WGYRE_PISCES_ST tracer.stat restartability passed : 14904
WORCA2_ICE_PISCES_ST run.stat restartability passed : 14904
WORCA2_ICE_PISCES_ST tracer.stat restartability passed : 14904
WORCA2_OFF_PISCES_ST tracer.stat restartability passed : 14904
WAMM12_ST run.stat restartability passed : 14904
WORCA2_SAS_ICE_ST run.stat restartability passed : 14904
WAGRIF_DEMO_ST run.stat restartability passed : 14904
WWED025_ST run.stat restartability passed : 14904
WISOMIP+_ST run.stat restartability passed : 14904
WOVERFLOW_ST run.stat restartability passed : 14904
WLOCK_EXCHANGE_ST run.stat restartability passed : 14904
WVORTEX_ST run.stat restartability passed : 14904
WICE_AGRIF_ST run.stat restartability passed : 14904
WSWG_ST run.stat restartability passed : 14904
!----repro----!
WGYRE_PISCES_ST run.stat reproducibility passed : 14904
WGYRE_PISCES_ST tracer.stat reproducibility passed : 14904
WORCA2_ICE_PISCES_ST run.stat reproducibility passed : 14904
WORCA2_ICE_PISCES_ST tracer.stat reproducibility passed : 14904
WORCA2_OFF_PISCES_ST tracer.stat reproducibility passed : 14904
WAMM12_ST run.stat reproducibility passed : 14904
WORCA2_SAS_ICE_ST run.stat reproducibility passed : 14904
WORCA2_ICE_OBS_ST run.stat reproducibility passed : 14904
WAGRIF_DEMO_ST run.stat reproducibility passed : 14904
WWED025_ST run.stat reproducibility passed : 14904
WISOMIP+_ST run.stat reproducibility passed : 14904
WVORTEX_ST run.stat reproducibility passed : 14904
WICE_AGRIF_ST run.stat reproducibility passed : 14904
WSWG_ST run.stat reproducibility passed : 14904
!----agrif check----!
ORCA2 AGRIF vs ORCA2 NOAGRIF run.stat unchanged - passed : 14904 14904
!----result comparison check----!
check result differences between :
VALID directory : /scratch/cnt0021/egi6035/pmathiot/NEMO/branches/ticket2669_isf_fluxes/NEMO_VALIDATION at rev 14904
and
REFERENCE directory : /scratch/cnt0021/egi6035/pmathiot/NEMO/trunk/NEMO_VALIDATION/ at rev 14886
WGYRE_PISCES_ST run.stat files are identical
WGYRE_PISCES_ST tracer.stat files are identical
WORCA2_ICE_PISCES_ST run.stat files are identical
WORCA2_ICE_PISCES_ST tracer.stat files are identical
WORCA2_OFF_PISCES_ST tracer.stat files are identical
WAMM12_ST run.stat files are identical
WORCA2_SAS_ICE_ST run.stat files are identical
WAGRIF_DEMO_ST run.stat files are identical
WWED025_ST run.stat files are DIFFERENT (results are different after 2 time steps)
WISOMIP+_ST run.stat files are identical
WVORTEX_ST run.stat files are identical
WICE_AGRIF_ST run.stat files are identical
WOVERFLOW_ST run.stat files are identical
WLOCK_EXCHANGE_ST run.stat files are identical
WSWG_ST run.stat files are identical
comment:15 Changed 3 years ago by mathiot
In 14926:
comment:16 Changed 3 years ago by mathiot
At the end, I did some simple tests with bg03, so I didn't add a STOP to prevent usage. I however did not change the doc on the quality (not so good) of bg03 itself.
comment:17 Changed 3 years ago by mathiot
At the end, I did some simple tests with bg03, so I didn't add a STOP to prevent usage. I however did not change the doc on the quality (not so good) of bg03 itself.
comment:18 Changed 3 years ago by mathiot
In 14988:
comment:19 Changed 3 years ago by mathiot
In 14989:
comment:20 Changed 3 years ago by mathiot
In 14994:
comment:21 Changed 3 years ago by mathiot
In 14995:
comment:22 Changed 3 years ago by mathiot
- Resolution set to fixed
- Status changed from assigned to closed
Erratum: I think issues are more dramatic than previously thought for closed cavities and parameterized melting. For bg03 and oasis, both latent and sensible fluxes have the wrong sign. One needs to double check but I think ocean is warming when isf melting instead of cooling
1) For specified melt, sensible flux has the wrong sign:
instead of
2) For bg03, both latent and sensible flux have the wrong sign (I think)
instead of
3) For oasis, both latent and sensible flux also have the wrong sign (I think)
instead of