Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1760 additions and 199 deletions
src/OCE/SBC/geo2ocean.F90
View file @ 4a9f4b18
......@@ -25,6 +25,7 @@ MODULE geo2ocean
IMPLICIT NONE
PRIVATE
PUBLIC repcmo ! called in sbccpl
PUBLIC rot_rep ! called in sbccpl, fldread, and cyclone
PUBLIC geo2oce ! called in sbccpl
PUBLIC oce2geo ! called in sbccpl
......@@ -50,6 +51,48 @@ MODULE geo2ocean
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1, &
px2 , py2 , kchoix )
!!----------------------------------------------------------------------
!! *** ROUTINE repcmo ***
!!
!! ** Purpose : Change vector componantes from a geographic grid to a
!! stretched coordinates grid.
!!
!! ** Method : Initialization of arrays at the first call.
!!
!! ** Action : - px2 : first componante (defined at u point)
!! - py2 : second componante (defined at v point)
!!----------------------------------------------------------------------
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxu1, pyu1 ! geographic vector componantes at u-point
REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxv1, pyv1 ! geographic vector componantes at v-point
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: px2 ! i-componante (defined at u-point)
REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: py2 ! j-componante (defined at v-point)
!!----------------------------------------------------------------------
INTEGER, INTENT( IN ) :: &
kchoix ! type of transformation
! = 1 change from geographic to model grid.
! =-1 change from model to geographic grid
!!----------------------------------------------------------------------
SELECT CASE (kchoix)
CASE ( 1)
! Change from geographic to stretched coordinate
! ----------------------------------------------
CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
CASE (-1)
! Change from stretched to geographic coordinate
! ----------------------------------------------
CALL rot_rep( pxu1, pyu1, 'U', 'ij->e',px2 )
CALL rot_rep( pxv1, pyv1, 'V', 'ij->n',py2 )
END SELECT
END SUBROUTINE repcmo
SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
!!----------------------------------------------------------------------
!! *** ROUTINE rot_rep ***
......
src/OCE/SBC/sbc_ice.F90
View file @ 4a9f4b18
......@@ -94,6 +94,7 @@ MODULE sbc_ice
! already defined in ice.F90 for SI3
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_i, h_s
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Sea ice fraction on categories at the last coupling point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tatm_ice !: air temperature [K]
#endif
......
src/OCE/SBC/sbc_oce.F90
View file @ 4a9f4b18
......@@ -131,6 +131,8 @@ MODULE sbc_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: atm_co2 !: atmospheric pCO2 [ppm]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: cloud_fra !: cloud cover (fraction of cloud in a gridcell) [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: greenland_icesheet_mask
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: antarctica_icesheet_mask
!!---------------------------------------------------------------------
!! ABL Vertical Domain size
......@@ -153,6 +155,33 @@ MODULE sbc_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3t_m !: mean (nn_fsbc time-step) sea surface layer thickness [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_m !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-]
!!----------------------------------------------------------------------
!! Surface scalars of total ice sheet mass for Greenland and Antarctica,
!! passed from atmosphere to be converted to dvol and hence a freshwater
!! flux by using old values. New values are saved in the dump, to become
!! old values next coupling timestep. Freshwater fluxes split between
!! sub iceshelf melting and iceberg calving, scalled to flux per second
!!----------------------------------------------------------------------
REAL(wp), PUBLIC :: greenland_icesheet_mass, greenland_icesheet_mass_rate_of_change, greenland_icesheet_timelapsed
REAL(wp), PUBLIC :: antarctica_icesheet_mass, antarctica_icesheet_mass_rate_of_change, antarctica_icesheet_timelapsed
! sbccpl namelist parameters associated with icesheet freshwater input code. Included here rather than in sbccpl.F90 to
! avoid circular dependencies.
INTEGER, PUBLIC :: nn_coupled_iceshelf_fluxes ! =0 : total freshwater input from iceberg calving and ice shelf basal melting
! taken from climatologies used (no action in coupling routines).
! =1 : use rate of change of mass of Greenland and Antarctic icesheets to set the
! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
! =2 : specify constant freshwater inputs in this namelist to set the combined
! magnitude of iceberg calving and iceshelf melting freshwater fluxes.
LOGICAL, PUBLIC :: ln_iceshelf_init_atmos ! If true force ocean to initialise iceshelf masses from atmospheric values rather
! than values in ocean restart (applicable if nn_coupled_iceshelf_fluxes=1).
REAL(wp), PUBLIC :: rn_greenland_total_fw_flux ! Constant total rate of freshwater input (kg/s) for Greenland (if nn_coupled_iceshelf_fluxes=2)
REAL(wp), PUBLIC :: rn_greenland_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
REAL(wp), PUBLIC :: rn_antarctica_total_fw_flux ! Constant total rate of freshwater input (kg/s) for Antarctica (if nn_coupled_iceshelf_fluxes=2)
REAL(wp), PUBLIC :: rn_antarctica_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
REAL(wp), PUBLIC :: rn_iceshelf_fluxes_tolerance ! Absolute tolerance for detecting differences in icesheet masses.
!!----------------------------------------------------------------------
!! Surface atmospheric fields
!!----------------------------------------------------------------------
......@@ -193,6 +222,9 @@ CONTAINS
& atm_co2(jpi,jpj) , tsk_m(jpi,jpj) , cloud_fra(jpi,jpj), &
& ssu_m (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) , &
& ssv_m (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) )
!
ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
!
......
src/OCE/SBC/sbcblk.F90
View file @ 4a9f4b18
......@@ -119,6 +119,7 @@ MODULE sbcblk
!
REAL(dp) :: rn_pfac ! multiplication factor for precipitation
REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation
REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress !
REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements
REAL(wp) :: rn_zu ! z(u) : height of wind measurements
!
......@@ -220,7 +221,7 @@ CONTAINS
INTEGER :: ipka ! number of levels in the atmospheric variable
NAMELIST/namsbc_blk/ ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, ln_ANDREAS, & ! bulk algorithm
& rn_zqt, rn_zu, nn_iter_algo, ln_skin_cs, ln_skin_wl, &
& rn_pfac, rn_efac, &
& rn_pfac, rn_efac, rn_vfac, &
& ln_crt_fbk, rn_stau_a, rn_stau_b, & ! current feedback
& ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tair_pot, &
& ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i, &
......@@ -414,6 +415,7 @@ CONTAINS
WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu
WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac
WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac
WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac
WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))'
WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk
IF(ln_crt_fbk) THEN
......@@ -655,9 +657,7 @@ CONTAINS
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zztmp ! local variable
REAL(wp) :: zstmax, zstau
#if defined key_cyclone
REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point
#endif
REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point
REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]
REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean
......@@ -683,9 +683,10 @@ CONTAINS
! ----------------------------------------------------------------------------- !
! ... components ( U10m - U_oce ) at T-point (unmasked)
#if defined key_cyclone
zwnd_i(:,:) = 0._wp
zwnd_j(:,:) = 0._wp
#if defined key_cyclone
CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
......@@ -695,9 +696,15 @@ CONTAINS
END_2D
#else
! ... scalar wind module at T-point (not masked)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
DO_2D( 0, 0, 0, 0 )
zwnd_i(ji,jj) = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )
zwnd_j(ji,jj) = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )
END_2D
CALL lbc_lnk( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &
& + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)
#endif
! ----------------------------------------------------------------------------- !
! I Solar FLUX !
......@@ -827,8 +834,13 @@ CONTAINS
ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
#else
ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
IF ( rn_vfac > 0._wp ) THEN
ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
ELSE
ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
ENDIF
#endif
ELSE
ztau_i(ji,jj) = 0._wp
......@@ -1017,7 +1029,8 @@ CONTAINS
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zootm_su ! sea-ice surface mean temperature
REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
REAL(wp) :: zwndi_t , zwndj_t
REAL(wp) :: zztmp0,zztmp1, zztmp2 ! temporary scalars
REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
!!---------------------------------------------------------------------
!
......@@ -1025,9 +1038,12 @@ CONTAINS
! Wind module relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
! C-grid ice dynamics : U & V-points (same as ocean)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
DO_2D( 0, 0, 0, 0 )
zwndi_t = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( puice(ji-1,jj ) + puice(ji,jj) ) )
zwndj_t = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pvice(ji ,jj-1) + pvice(ji,jj) ) )
wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t )* tmask(ji,jj,1)
END_2D
CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. )
!
! potential sea-ice surface temperature [K]
zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
......@@ -1074,10 +1090,11 @@ CONTAINS
! Wind stress relative to nonmoving ice ( U10m ) !
! ---------------------------------------------------- !
! supress moving ice in wind stress computation as we don't know how to do it properly...
zztmp0 = rn_vfac * 0.5_wp
DO_2D( 0, 1, 0, 1 ) ! at T point
zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
putaui(ji,jj) = zztmp1 * ( pwndi(ji,jj) - zztmp0 * ( puice(ji-1,jj ) + puice(ji,jj) ) )
pvtaui(ji,jj) = zztmp1 * ( pwndj(ji,jj) - zztmp0 * ( pvice(ji ,jj-1) + pvice(ji,jj) ) )
END_2D
!#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!!
......
src/OCE/SBC/sbccpl.F90
View file @ 4a9f4b18
......@@ -33,10 +33,18 @@ MODULE sbccpl
#endif
USE cpl_oasis3 ! OASIS3 coupling
USE geo2ocean !
USE oce , ONLY : ts, uu, vv, ssh, fraqsr_1lev
USE oce
#if defined key_medusa
USE oce , ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl, &
PCO2a_in_cpl, Dust_in_cpl
#endif
USE ocealb !
USE eosbn2 !
USE sbcrnf , ONLY : l_rnfcpl
USE cpl_rnf_1d, ONLY: nn_cpl_river, cpl_rnf_1d_init, cpl_rnf_1d_to_2d ! Variables used in 1D river outflow
#if defined key_medusa
USE par_trc , ONLY : ln_medusa
#endif
#if defined key_cice
USE ice_domain_size, only: ncat
#endif
......@@ -48,6 +56,7 @@ MODULE sbccpl
USE iom ! NetCDF library
USE lib_mpp ! distribued memory computing library
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE timing
#if defined key_oasis3
USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut
......@@ -129,8 +138,17 @@ MODULE sbccpl
INTEGER, PARAMETER :: jpr_icb = 61
INTEGER, PARAMETER :: jpr_ts_ice = 62 ! Sea ice surface temp
!!INTEGER, PARAMETER :: jpr_qtrice = 63 ! Transmitted solar thru sea-ice
INTEGER, PARAMETER :: jprcv = 62 ! total number of fields received
INTEGER, PARAMETER :: jpr_grnm = 63 ! Greenland ice mass
INTEGER, PARAMETER :: jpr_antm = 64 ! Antarctic ice mass
INTEGER, PARAMETER :: jpr_rnf_1d = 65 ! 1D river runoff
INTEGER, PARAMETER :: jpr_qtr = 66 ! Transmitted solar
#if defined key_medusa
INTEGER, PARAMETER :: jpr_atm_pco2 = 67 ! Incoming atm pCO2 flux
INTEGER, PARAMETER :: jpr_atm_dust = 68 ! Incoming atm aggregate dust
INTEGER, PARAMETER :: jprcv = 69 ! total number of fields received
#else
INTEGER, PARAMETER :: jprcv = 66 ! total number of fields received
#endif
INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere
INTEGER, PARAMETER :: jps_toce = 2 ! ocean temperature
......@@ -170,17 +188,24 @@ MODULE sbccpl
INTEGER, PARAMETER :: jps_kice = 36 ! sea ice effective conductivity
INTEGER, PARAMETER :: jps_sstfrz = 37 ! sea surface freezing temperature
INTEGER, PARAMETER :: jps_ttilyr = 38 ! sea ice top layer temp
INTEGER, PARAMETER :: jpsnd = 38 ! total number of fields sent
#if ! defined key_oasis3
! Dummy variables to enable compilation when oasis3 is not being used
INTEGER :: OASIS_Sent = -1
INTEGER :: OASIS_SentOut = -1
INTEGER :: OASIS_ToRest = -1
INTEGER :: OASIS_ToRestOut = -1
#if defined key_medusa
INTEGER, PARAMETER :: jps_bio_co2 = 39 ! MEDUSA air-sea CO2 flux
INTEGER, PARAMETER :: jps_bio_dms = 40 ! MEDUSA DMS surface concentration
INTEGER, PARAMETER :: jps_bio_chloro = 41 ! MEDUSA chlorophyll surface concentration
INTEGER, PARAMETER :: jpsnd = 41 ! total number of fields sent
#else
INTEGER, PARAMETER :: jpsnd = 38 ! total number of fields sent
#endif
#if ! defined key_oasis3
! Dummy variables to enable compilation when oasis3 is not being used
INTEGER :: OASIS_Sent = -1
INTEGER :: OASIS_SentOut = -1
INTEGER :: OASIS_ToRest = -1
INTEGER :: OASIS_ToRestOut = -1
#endif
! !!** namelist namsbc_cpl **
TYPE :: FLD_C !
CHARACTER(len = 32) :: cldes ! desciption of the coupling strategy
......@@ -193,20 +218,34 @@ MODULE sbccpl
TYPE(FLD_C) :: sn_snd_temp , sn_snd_alb , sn_snd_thick, sn_snd_crt , sn_snd_co2, &
& sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
! ! Received from the atmosphere
TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, &
#if defined key_medusa
TYPE(FLD_C) :: sn_snd_bio_co2, sn_snd_bio_dms, sn_snd_bio_chloro
#endif
TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr, &
& sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf, sn_rcv_ts_ice
TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
! ! Send to waves
TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev
! ! Received from waves
TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf, &
sn_rcv_grnm, sn_rcv_antm
#if defined key_medusa
TYPE(FLD_C) :: sn_rcv_atm_pco2, sn_rcv_atm_dust
#endif
! Send to waves
TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev
TYPE(FLD_C) :: sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, &
& sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
! Received from waves
TYPE(FLD_C) :: sn_rcv_tauwoc, sn_rcv_wfreq
! Transmitted solar
TYPE(FLD_C) :: sn_rcv_qtr
! ! Other namelist parameters
!! TYPE(FLD_C) :: sn_rcv_qtrice
INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
LOGICAL :: ln_usecplmask ! use a coupling mask file to merge data received from several models
! -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
LOGICAL :: ln_scale_ice_flux ! use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)
LOGICAL :: ln_couple_ocean_evap ! Do we couple total (ocean+sea ice) evaporation (FALSE)
! or ocean only evaporation (TRUE)
TYPE :: DYNARR
REAL(wp), POINTER, DIMENSION(:,:,:) :: z3
......@@ -215,9 +254,7 @@ MODULE sbccpl
TYPE( DYNARR ), SAVE, DIMENSION(jprcv) :: frcv ! all fields recieved from the atmosphere
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: alb_oce_mix ! ocean albedo sent to atmosphere (mix clear/overcast sky)
#if defined key_si3 || defined key_cice
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Ice fractional area at last coupling time
#endif
INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) :: nrcvinfo ! OASIS info argument
......@@ -275,22 +312,38 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: zacs, zaos
!!
NAMELIST/namsbc_cpl/ nn_cplmodel , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux, &
& ln_couple_ocean_evap, &
& sn_snd_temp , sn_snd_alb , sn_snd_thick, sn_snd_crt , sn_snd_co2 , &
& sn_snd_ttilyr, sn_snd_cond , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1, &
& sn_snd_ifrac , sn_snd_crtw , sn_snd_wlev , sn_rcv_hsig , sn_rcv_phioc , &
& sn_rcv_w10m , sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr , &
& sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum , sn_rcv_wstrf , &
& sn_rcv_charn , sn_rcv_taw , sn_rcv_bhd , sn_rcv_tusd , sn_rcv_tvsd, &
sn_rcv_qtr , &
& sn_rcv_wdrag , sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , &
& sn_rcv_iceflx, sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf , sn_rcv_ts_ice, & !!, sn_rcv_qtrice
& sn_rcv_mslp
& sn_rcv_mslp, &
sn_rcv_grnm , sn_rcv_antm , &
& nn_coupled_iceshelf_fluxes , ln_iceshelf_init_atmos , &
& rn_greenland_total_fw_flux , rn_greenland_calving_fraction , &
& rn_antarctica_total_fw_flux , rn_antarctica_calving_fraction , &
#if defined key_medusa
& rn_iceshelf_fluxes_tolerance, &
! Add MEDUSA related fields to namelist
sn_snd_bio_co2 , sn_snd_bio_dms, sn_snd_bio_chloro, &
& sn_rcv_atm_pco2, sn_rcv_atm_dust
#else
& rn_iceshelf_fluxes_tolerance
#endif
!!---------------------------------------------------------------------
!
! ================================ !
! Namelist informations !
! ================================ !
!
!
IF (ln_timing) CALL timing_start('sbc_cpl_init')
READ ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist' )
!
......@@ -307,6 +360,7 @@ CONTAINS
WRITE(numout,*)' nn_cplmodel = ', nn_cplmodel
WRITE(numout,*)' ln_usecplmask = ', ln_usecplmask
WRITE(numout,*)' ln_scale_ice_flux = ', ln_scale_ice_flux
WRITE(numout,*)' ln_couple_ocean_evap = ', ln_couple_ocean_evap
WRITE(numout,*)' nn_cats_cpl = ', nn_cats_cpl
WRITE(numout,*)' received fields (mutiple ice categogies)'
WRITE(numout,*)' 10m wind module = ', TRIM(sn_rcv_w10m%cldes ), ' (', TRIM(sn_rcv_w10m%clcat ), ')'
......@@ -321,10 +375,12 @@ CONTAINS
WRITE(numout,*)' freshwater budget = ', TRIM(sn_rcv_emp%cldes ), ' (', TRIM(sn_rcv_emp%clcat ), ')'
WRITE(numout,*)' runoffs = ', TRIM(sn_rcv_rnf%cldes ), ' (', TRIM(sn_rcv_rnf%clcat ), ')'
WRITE(numout,*)' calving = ', TRIM(sn_rcv_cal%cldes ), ' (', TRIM(sn_rcv_cal%clcat ), ')'
WRITE(numout,*)' Greenland ice mass = ', TRIM(sn_rcv_grnm%cldes ), ' (', TRIM(sn_rcv_grnm%clcat ), ')'
WRITE(numout,*)' Antarctica ice mass = ', TRIM(sn_rcv_antm%cldes ), ' (', TRIM(sn_rcv_antm%clcat ), ')'
WRITE(numout,*)' iceberg = ', TRIM(sn_rcv_icb%cldes ), ' (', TRIM(sn_rcv_icb%clcat ), ')'
WRITE(numout,*)' ice shelf = ', TRIM(sn_rcv_isf%cldes ), ' (', TRIM(sn_rcv_isf%clcat ), ')'
WRITE(numout,*)' sea ice heat fluxes = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
!! WRITE(numout,*)' transmitted solar thru sea-ice = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%clcat), ')'
WRITE(numout,*)' transmitted solar = ', TRIM(sn_rcv_qtr%cldes ), ' (', TRIM(sn_rcv_qtr%clcat ), ')'
WRITE(numout,*)' atm co2 = ', TRIM(sn_rcv_co2%cldes ), ' (', TRIM(sn_rcv_co2%clcat ), ')'
WRITE(numout,*)' Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'
WRITE(numout,*)' surface waves:'
......@@ -347,6 +403,11 @@ CONTAINS
WRITE(numout,*)' - referential = ', sn_snd_crt%clvref
WRITE(numout,*)' - orientation = ', sn_snd_crt%clvor
WRITE(numout,*)' - mesh = ', sn_snd_crt%clvgrd
#if defined key_medusa
WRITE(numout,*)' bio co2 flux = ', TRIM(sn_snd_bio_co2%cldes),' (', TRIM(sn_snd_bio_co2%clcat), ')'
WRITE(numout,*)' bio dms flux = ', TRIM(sn_snd_bio_dms%cldes),' (', TRIM(sn_snd_bio_dms%clcat), ')'
WRITE(numout,*)' bio dms chlorophyll = ', TRIM(sn_snd_bio_chloro%cldes), ' (', TRIM(sn_snd_bio_chloro%clcat), ')'
#endif
WRITE(numout,*)' oce co2 flux = ', TRIM(sn_snd_co2%cldes ), ' (', TRIM(sn_snd_co2%clcat ), ')'
WRITE(numout,*)' ice effective conductivity = ', TRIM(sn_snd_cond%cldes ), ' (', TRIM(sn_snd_cond%clcat ), ')'
WRITE(numout,*)' meltponds fraction and depth = ', TRIM(sn_snd_mpnd%cldes ), ' (', TRIM(sn_snd_mpnd%clcat ), ')'
......@@ -357,6 +418,17 @@ CONTAINS
WRITE(numout,*)' - referential = ', sn_snd_crtw%clvref
WRITE(numout,*)' - orientation = ', sn_snd_crtw%clvor
WRITE(numout,*)' - mesh = ', sn_snd_crtw%clvgrd
#if defined key_medusa
WRITE(numout,*)' atm pco2 = ', TRIM(sn_rcv_atm_pco2%cldes),'(', TRIM(sn_rcv_atm_pco2%clcat ), ')'
WRITE(numout,*)' atm dust = ', TRIM(sn_rcv_atm_dust%cldes),'(', TRIM(sn_rcv_atm_dust%clcat),')'
#endif
WRITE(numout,*)' nn_coupled_iceshelf_fluxes = ', nn_coupled_iceshelf_fluxes
WRITE(numout,*)' ln_iceshelf_init_atmos = ', ln_iceshelf_init_atmos
WRITE(numout,*)' rn_greenland_total_fw_flux = ', rn_greenland_total_fw_flux
WRITE(numout,*)' rn_antarctica_total_fw_flux = ', rn_antarctica_total_fw_flux
WRITE(numout,*)' rn_greenland_calving_fraction = ', rn_greenland_calving_fraction
WRITE(numout,*)' rn_antarctica_calving_fraction = ', rn_antarctica_calving_fraction
WRITE(numout,*)' rn_iceshelf_fluxes_tolerance = ', rn_iceshelf_fluxes_tolerance
ENDIF
IF( lwp .AND. ln_wave) THEN ! control print
WRITE(numout,*)' surface waves:'
......@@ -391,7 +463,11 @@ CONTAINS
! define the north fold type of lbc (srcv(:)%nsgn)
! default definitions of srcv
srcv(:)%laction = .FALSE. ; srcv(:)%clgrid = 'T' ; srcv(:)%nsgn = 1. ; srcv(:)%nct = 1
srcv(:)%laction = .FALSE.
srcv(:)%clgrid = 'T'
srcv(:)%nsgn = 1.
srcv(:)%nct = 1
srcv(:)%dimensions = 2
! ! ------------------------- !
! ! ice and ocean wind stress !
......@@ -445,7 +521,11 @@ CONTAINS
srcv(jpr_otx1:jpr_otz1)%clgrid = 'U' ! oce components given at U-point
srcv(jpr_otx2:jpr_otz2)%clgrid = 'V' ! and V-point
srcv(jpr_itx1:jpr_itz1)%clgrid = 'F' ! ice components given at F-point
srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
!srcv(jpr_otx1:jpr_otz2)%laction = .TRUE. ! receive oce components on grid 1 & 2
! Currently needed for HadGEM3 - but shouldn't affect anyone else for the moment
srcv(jpr_otx1)%laction = .TRUE.
srcv(jpr_oty1)%laction = .TRUE.
!
srcv(jpr_itx1:jpr_itz1)%laction = .TRUE. ! receive ice components on grid 1 only
CASE( 'T,I' )
srcv(jpr_otx1:jpr_itz2)%clgrid = 'T' ! oce and ice components given at T-point
......@@ -510,15 +590,34 @@ CONTAINS
! ! Runoffs & Calving !
! ! ------------------------- !
srcv(jpr_rnf )%clname = 'O_Runoff'
IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
srcv(jpr_rnf)%laction = .TRUE.
srcv(jpr_rnf_1d )%clname = 'ORunff1D'
IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' .OR. TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN
IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) srcv(jpr_rnf)%laction = .TRUE.
IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN
srcv(jpr_rnf_1d)%laction = .TRUE.
srcv(jpr_rnf_1d)%dimensions = 1 ! 1D field passed through coupler
END IF
l_rnfcpl = .TRUE. ! -> no need to read runoffs in sbcrnf
ln_rnf = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' runoffs received from oasis -> force ln_rnf = ', ln_rnf
ENDIF
!
srcv(jpr_cal)%clname = 'OCalving' ; IF( TRIM( sn_rcv_cal%cldes) == 'coupled' ) srcv(jpr_cal)%laction = .TRUE.
srcv(jpr_cal )%clname = 'OCalving'
IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' ) srcv(jpr_cal)%laction = .TRUE.
srcv(jpr_grnm )%clname = 'OGrnmass'
IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' ) THEN
srcv(jpr_grnm)%laction = .TRUE.
srcv(jpr_grnm)%dimensions = 0 ! Scalar field
ENDIF
srcv(jpr_antm )%clname = 'OAntmass'
IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' ) THEN
srcv(jpr_antm)%laction = .TRUE.
srcv(jpr_antm)%dimensions = 0 ! Scalar field
ENDIF
srcv(jpr_isf)%clname = 'OIcshelf' ; IF( TRIM( sn_rcv_isf%cldes) == 'coupled' ) srcv(jpr_isf)%laction = .TRUE.
srcv(jpr_icb)%clname = 'OIceberg' ; IF( TRIM( sn_rcv_icb%cldes) == 'coupled' ) srcv(jpr_icb)%laction = .TRUE.
......@@ -594,6 +693,22 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' Atmospheric pco2 received from oasis '
IF(lwp) WRITE(numout,*)
ENDIF
#if defined key_medusa
! ! --------------------------------------- !
! ! Incoming CO2 and DUST fluxes for MEDUSA !
! ! --------------------------------------- !
srcv(jpr_atm_pco2)%clname = 'OATMPCO2'
IF (TRIM(sn_rcv_atm_pco2%cldes) == 'medusa') THEN
srcv(jpr_atm_pco2)%laction = .TRUE.
END IF
srcv(jpr_atm_dust)%clname = 'OATMDUST'
IF (TRIM(sn_rcv_atm_dust%cldes) == 'medusa') THEN
srcv(jpr_atm_dust)%laction = .TRUE.
END IF
#endif
!
! ! ------------------------- !
! ! Mean Sea Level Pressure !
......@@ -625,6 +740,21 @@ CONTAINS
!! ENDIF
!! srcv(jpr_qtrice)%laction = .TRUE.
!! ENDIF
! ! ------------------------- !
! ! transmitted solar !
! ! ------------------------- !
srcv(jpr_qtr )%clname = 'OQtr'
IF( TRIM(sn_rcv_qtr%cldes) == 'coupled' ) THEN
IF ( TRIM( sn_rcv_qtr%clcat ) == 'yes' ) THEN
srcv(jpr_qtr)%nct = nn_cats_cpl
ELSE
CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtr%clcat should always be set to yes currently' )
ENDIF
srcv(jpr_qtr)%laction = .TRUE.
ENDIF
! ! ------------------------- !
! ! ice skin temperature !
! ! ------------------------- !
......@@ -799,24 +929,6 @@ CONTAINS
ENDIF
ENDIF
! =================================================== !
! Allocate all parts of frcv used for received fields !
! =================================================== !
DO jn = 1, jprcv
IF( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
END DO
! Allocate taum part of frcv which is used even when not received as coupling field
IF( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
! Allocate w10m part of frcv which is used even when not received as coupling field
IF( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
IF( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
IF( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
IF( k_ice /= 0 ) THEN
IF( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
IF( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
ENDIF
! ================================ !
! Define the send interface !
......@@ -826,8 +938,12 @@ CONTAINS
! define the north fold type of lbc (ssnd(:)%nsgn)
! default definitions of nsnd
ssnd(:)%laction = .FALSE. ; ssnd(:)%clgrid = 'T' ; ssnd(:)%nsgn = 1. ; ssnd(:)%nct = 1
ssnd(:)%laction = .FALSE.
ssnd(:)%clgrid = 'T'
ssnd(:)%nsgn = 1.
ssnd(:)%nct = 1
ssnd(:)%dimensions = 2
! ! ------------------------- !
! ! Surface temperature !
! ! ------------------------- !
......@@ -978,6 +1094,25 @@ CONTAINS
! ! ------------------------- !
ssnd(jps_co2)%clname = 'O_CO2FLX' ; IF( TRIM(sn_snd_co2%cldes) == 'coupled' ) ssnd(jps_co2 )%laction = .TRUE.
!
#if defined key_medusa
! ! ------------------------- !
! ! MEDUSA output fields !
! ! ------------------------- !
! Surface dimethyl sulphide from Medusa
ssnd(jps_bio_dms)%clname = 'OBioDMS'
IF( TRIM(sn_snd_bio_dms%cldes) == 'medusa' ) ssnd(jps_bio_dms )%laction = .TRUE.
! Surface CO2 flux from Medusa
ssnd(jps_bio_co2)%clname = 'OBioCO2'
IF( TRIM(sn_snd_bio_co2%cldes) == 'medusa' ) ssnd(jps_bio_co2 )%laction = .TRUE.
! Surface chlorophyll from Medusa
ssnd(jps_bio_chloro)%clname = 'OBioChlo'
IF( TRIM(sn_snd_bio_chloro%cldes) == 'medusa' ) ssnd(jps_bio_chloro )%laction = .TRUE.
#endif
! ! ------------------------- !
! ! Sea surface freezing temp !
! ! ------------------------- !
......@@ -1105,10 +1240,52 @@ CONTAINS
ENDIF
ENDIF
! Initialise 1D river outflow scheme
nn_cpl_river = 1
IF ( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) CALL cpl_rnf_1d_init ! Coupled runoff using 1D array
! =================================================== !
! Allocate all parts of frcv used for received fields !
! =================================================== !
DO jn = 1, jprcv
IF ( srcv(jn)%laction ) THEN
SELECT CASE( srcv(jn)%dimensions )
!
CASE( 0 ) ! Scalar field
ALLOCATE( frcv(jn)%z3(1,1,1) )
CASE( 1 ) ! 1D field
ALLOCATE( frcv(jn)%z3(nn_cpl_river,1,1) )
CASE DEFAULT ! 2D (or pseudo 3D) field.
ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
END SELECT
END IF
END DO
! Allocate taum part of frcv which is used even when not received as coupling field
IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
! Allocate w10m part of frcv which is used even when not received as coupling field
IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
IF( k_ice /= 0 ) THEN
IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
END IF
!
! ================================ !
! initialisation of the coupler !
! ================================ !
! There's no point initialising the coupler if we've accumulated any errors in
! coupling field definitions or settings.
IF (nstop > 0) CALL ctl_stop( 'STOP', 'sbc_cpl_init: Errors encountered in coupled field definitions' )
CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
IF(ln_usecplmask) THEN
......@@ -1122,6 +1299,26 @@ CONTAINS
ENDIF
xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
!
IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something
! more complicated could be done if required.
greenland_icesheet_mask = 0.0
WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0
antarctica_icesheet_mask = 0.0
WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0
IF( .not. ln_rstart ) THEN
greenland_icesheet_mass = 0.0
greenland_icesheet_mass_rate_of_change = 0.0
greenland_icesheet_timelapsed = 0.0
antarctica_icesheet_mass = 0.0
antarctica_icesheet_mass_rate_of_change = 0.0
antarctica_icesheet_timelapsed = 0.0
ENDIF
ENDIF
!
IF (ln_timing) CALL timing_stop('sbc_cpl_init')
!
END SUBROUTINE sbc_cpl_init
......@@ -1181,15 +1378,26 @@ CONTAINS
LOGICAL :: llnewtx, llnewtau ! update wind stress components and module??
INTEGER :: ji, jj, jn ! dummy loop indices
INTEGER :: isec ! number of seconds since nit000 (assuming rdt did not change since nit000)
INTEGER :: ikchoix
REAL(wp), DIMENSION(jpi,jpj) :: ztx2, zty2
REAL(wp) :: zcumulneg, zcumulpos ! temporary scalars
REAL(wp) :: zcoef ! temporary scalar
LOGICAL :: ll_wrtstp ! write diagnostics?
REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
REAL(wp) :: zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in
REAL(wp) :: zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b
REAL(wp) :: zmask_sum, zepsilon
REAL(wp) :: zzx, zzy ! temporary variables
REAL(wp) :: r1_grau ! = 1.e0 / (grav * rho0)
REAL(wp), DIMENSION(jpi,jpj) :: ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
!!----------------------------------------------------------------------
!
!
IF (ln_timing) CALL timing_start('sbc_cpl_rcv')
!
ll_wrtstp = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
!
IF( kt == nit000 ) THEN
! cannot be done in the init phase when we use agrif as cpl_freq requires that oasis_enddef is done
ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
......@@ -1209,7 +1417,16 @@ CONTAINS
! ! ======================================================= !
isec = ( kt - nit000 ) * NINT( rn_Dt ) ! date of exchanges
DO jn = 1, jprcv ! received fields sent by the atmosphere
IF( srcv(jn)%laction ) CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
IF( srcv(jn)%laction ) THEN
IF ( srcv(jn)%dimensions <= 1 ) THEN
CALL cpl_rcv_1d( jn, isec, frcv(jn)%z3, SIZE(frcv(jn)%z3), nrcvinfo(jn) )
ELSE
CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
END IF
END IF
END DO
! ! ========================= !
......@@ -1238,14 +1455,42 @@ CONTAINS
!
IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN ! 2 components oriented along the local grid
! ! (geographical to local grid -> rotate the components)
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )
IF( srcv(jpr_otx2)%laction ) THEN
CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )
IF( srcv(jpr_otx1)%clgrid == 'U' .AND. (.NOT. srcv(jpr_otx2)%laction) ) THEN
! Temporary code for HadGEM3 - will be removed eventually.
! Only applies when we have only taux on U grid and tauy on V grid
!RSRH these MUST be initialised because the halos are not explicitly set
! but they're passed to repcmo and used directly in calculations, so if
! they point at junk in memory then bad things will happen!
! (You can prove this by running with preset NaNs).
ztx(:,:)=0.0
zty(:,:)=0.0
DO_2D( 0, 0, 0, 0 )
ztx(ji,jj)=0.25*vmask(ji,jj,1) &
*(frcv(jpr_otx1)%z3(ji,jj,1)+frcv(jpr_otx1)%z3(ji-1,jj,1) &
+frcv(jpr_otx1)%z3(ji,jj+1,1)+frcv(jpr_otx1)%z3(ji-1,jj+1,1))
zty(ji,jj)=0.25*umask(ji,jj,1) &
*(frcv(jpr_oty1)%z3(ji,jj,1)+frcv(jpr_oty1)%z3(ji+1,jj,1) &
+frcv(jpr_oty1)%z3(ji,jj-1,1)+frcv(jpr_oty1)%z3(ji+1,jj-1,1))
END_2D
ikchoix = 1
CALL repcmo (frcv(jpr_otx1)%z3(:,:,1),zty,ztx,frcv(jpr_oty1)%z3(:,:,1),ztx2,zty2,ikchoix)
CALL lbc_lnk ('jpr_otx1', ztx2,'U', -1. )
CALL lbc_lnk ('jpr_oty1', zty2,'V', -1. )
frcv(jpr_otx1)%z3(:,:,1)=ztx2(:,:)
frcv(jpr_oty1)%z3(:,:,1)=zty2(:,:)
ELSE
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )
frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st component on the 1st grid
IF( srcv(jpr_otx2)%laction ) THEN
CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )
ELSE
CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
ENDIF
frcv(jpr_oty1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd component on the 2nd grid
ENDIF
frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:) ! overwrite 1st component on the 1st grid
frcv(jpr_oty1)%z3(:,:,1) = zty(:,:) ! overwrite 2nd component on the 2nd grid
ENDIF
!
IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
......@@ -1335,6 +1580,12 @@ CONTAINS
!
ENDIF
#if defined key_medusa
IF (ln_medusa) THEN
IF( srcv(jpr_atm_pco2)%laction) PCO2a_in_cpl(:,:) = frcv(jpr_atm_pco2)%z3(:,:,1)
IF( srcv(jpr_atm_dust)%laction) Dust_in_cpl(:,:) = frcv(jpr_atm_dust)%z3(:,:,1)
ENDIF
#endif
! ! ================== !
! ! atmosph. CO2 (ppm) !
! ! ================== !
......@@ -1539,8 +1790,10 @@ CONTAINS
!
IF( srcv(jpr_icb)%laction ) zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus ! remove heat content associated to iceberg melting
!
IF( ln_mixcpl ) THEN ; qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
ELSE ; qns(:,:) = zqns(:,:)
IF( ln_mixcpl ) THEN
qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
ELSE
qns(:,:) = zqns(:,:)
ENDIF
! ! solar flux over the ocean (qsr)
......@@ -1559,7 +1812,78 @@ CONTAINS
IF( srcv(jpr_fice )%laction ) fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
!
ENDIF
zepsilon = rn_iceshelf_fluxes_tolerance
IF( srcv(jpr_grnm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
! This is a zero dimensional, single value field.
zgreenland_icesheet_mass_in = frcv(jpr_grnm)%z3(1,1,1)
greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt
IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
zgreenland_icesheet_mass_b = zgreenland_icesheet_mass_in
greenland_icesheet_mass = zgreenland_icesheet_mass_in
ENDIF
IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN
zgreenland_icesheet_mass_b = greenland_icesheet_mass
! Only update the mass if it has increased.
IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN
greenland_icesheet_mass = zgreenland_icesheet_mass_in
ENDIF
IF( zgreenland_icesheet_mass_b /= 0.0 ) &
& greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed
greenland_icesheet_timelapsed = 0.0_wp
ENDIF
IF(lwp .AND. ll_wrtstp) THEN
WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
WRITE(numout,*) 'Greenland icesheet mass (kg) used is ', greenland_icesheet_mass
WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
ENDIF
ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
greenland_icesheet_mass_rate_of_change = rn_greenland_total_fw_flux
ENDIF
! ! land ice masses : Antarctica
IF( srcv(jpr_antm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
! This is a zero dimensional, single value field.
zantarctica_icesheet_mass_in = frcv(jpr_antm)%z3(1,1,1)
antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt
IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
zantarctica_icesheet_mass_b = zantarctica_icesheet_mass_in
antarctica_icesheet_mass = zantarctica_icesheet_mass_in
ENDIF
IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN
zantarctica_icesheet_mass_b = antarctica_icesheet_mass
! Only update the mass if it has increased.
IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN
antarctica_icesheet_mass = zantarctica_icesheet_mass_in
END IF
IF( zantarctica_icesheet_mass_b /= 0.0 ) &
& antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed
antarctica_icesheet_timelapsed = 0.0_wp
ENDIF
IF(lwp .AND. ll_wrtstp) THEN
WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
WRITE(numout,*) 'Antarctica icesheet mass (kg) used is ', antarctica_icesheet_mass
WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
ENDIF
ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
ENDIF
!
IF (ln_timing) CALL timing_stop('sbc_cpl_rcv')
END SUBROUTINE sbc_cpl_rcv
......@@ -1759,17 +2083,32 @@ CONTAINS
!!----------------------------------------------------------------------
!
#if defined key_si3 || defined key_cice
!
IF( kt == nit000 ) THEN
! allocate ice fractions from last coupling time here and not in sbc_cpl_init because of jpl
IF( .NOT.ALLOCATED(a_i_last_couple) ) ALLOCATE( a_i_last_couple(jpi,jpj,jpl) )
! initialize to a_i for the 1st time step
a_i_last_couple(:,:,:) = a_i(:,:,:)
ENDIF
!
IF( ln_mixcpl ) zmsk(:,:) = 1. - xcplmask(:,:,0)
ziceld(:,:) = 1._wp - picefr(:,:)
zcptn (:,:) = rcp * sst_m(:,:)
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) ) ! ice topmelt and botmelt !
! ! ========================= !
CASE ('coupled')
IF (ln_scale_ice_flux) THEN
WHERE( a_i(:,:,:) > 1.e-10_wp )
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
ELSEWHERE
qml_ice(:,:,:) = 0.0_wp
qcn_ice(:,:,:) = 0.0_wp
END WHERE
ELSE
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
ENDIF
END SELECT
!
! ! ========================= !
! ! freshwater budget ! (emp_tot)
......@@ -1783,7 +2122,9 @@ CONTAINS
CASE( 'conservative' ) ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1) ! May need to ensure positive here
ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:) ! May need to ensure positive here
zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
IF (.not. ln_couple_ocean_evap ) THEN
zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
END IF
CASE( 'oce and ice' ) ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * picefr(:,:)
......@@ -1791,7 +2132,9 @@ CONTAINS
ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
CASE( 'none' ) ! Not available as for now: needs additional coding below when computing zevap_oce
! ! since fields received are not defined with none option
CALL ctl_stop('STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl')
CALL ctl_stop( 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
CASE default ! Default
CALL ctl_stop( 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
END SELECT
! --- evaporation over ice (kg/m2/s) --- !
......@@ -1837,19 +2180,30 @@ CONTAINS
! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( picefr(:,:) - zsnw(:,:) ) ! emp_ice = A * sublimation - zsnw * sprecip
zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) ! emp_oce = emp_tot - emp_ice
! --- evaporation over ocean (used later for qemp) --- !
zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
IF ( ln_couple_ocean_evap ) THEN
zemp_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_rain)%z3(:,:,1) & !Ocean evap minus rain (as all rain goes straight to ocean in GC5)
& - zsprecip(:,:) * ( 1._wp - zsnw(:,:) ) !subtract snow in leads after correction for blowing snow
zemp_tot(:,:) = zemp_oce(:,:) + zemp_ice(:,:)
zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1)
ELSE
zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) ! emp_oce = emp_tot - emp_ice
! --- evaporation over ocean (used later for qemp) --- !
zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
END IF
! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
! therefore, sublimation is not redistributed over the ice categories when no subgrid scale fluxes are provided by atm.
zdevap_ice(:,:) = 0._wp
! --- Continental fluxes --- !
IF( srcv(jpr_rnf)%laction ) THEN ! runoffs (included in emp later on)
IF( srcv(jpr_rnf)%laction ) THEN ! 2D runoffs (included in emp later on)
rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
ENDIF
IF( srcv(jpr_rnf_1d)%laction ) THEN ! 1D runoff
CALL cpl_rnf_1d_to_2d(frcv(jpr_rnf_1d)%z3(:,:,:))
ENDIF
IF( srcv(jpr_cal)%laction ) THEN ! calving (put in emp_tot and emp_oce)
zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
......@@ -1929,23 +2283,7 @@ CONTAINS
!! IF( srcv(jpr_rnf)%laction ) CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1) ) ! runoff
!! IF( srcv(jpr_isf)%laction ) CALL iom_put( 'iceshelf_cea', fwfisf(:,:) * tmask(:,:,1) ) ! iceshelf
!
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) ) ! ice topmelt and botmelt !
! ! ========================= !
CASE ('coupled')
IF (ln_scale_ice_flux) THEN
WHERE( a_i(:,:,:) > 1.e-10_wp )
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
ELSEWHERE
qml_ice(:,:,:) = 0.0_wp
qcn_ice(:,:,:) = 0.0_wp
END WHERE
ELSE
qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
ENDIF
END SELECT
!
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_qns%cldes ) ) ! non solar heat fluxes ! (qns)
......@@ -2009,8 +2347,10 @@ CONTAINS
! --- non solar flux over ocean --- !
! note: ziceld cannot be = 0 since we limit the ice concentration to amax
zqns_oce = 0._wp
WHERE( ziceld /= 0._wp ) zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
WHERE( ziceld /= 0._wp ) zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
! Heat content per unit mass of snow (J/kg)
WHERE( SUM( a_i, dim=3 ) > 1.e-10 ) ; zcptsnw(:,:) = rcpi * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 )
ELSEWHERE ; zcptsnw(:,:) = zcptn(:,:)
......@@ -2023,7 +2363,12 @@ CONTAINS
! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
DO jl = 1, jpl
zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account
! zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account
! How much enthalpy is stored in sublimating snow and ice. At this stage we don't know if it is snow or ice that is
! sublimating so we will use the combined snow and ice layer temperature t1_ice.
zqevap_ice(:,:,jl) = -zevap_ice(:,:,jl) * ( ( rt0 - t1_ice(:,:,jl) ) * rcpi + rLfus )
END DO
! --- heat flux associated with emp (W/m2) --- !
......@@ -2041,6 +2386,7 @@ CONTAINS
IF( ln_mixcpl ) THEN
qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
DO jl=1,jpl
qns_ice (:,:,jl) = qns_ice (:,:,jl) * xcplmask(:,:,0) + zqns_ice (:,:,jl)* zmsk(:,:)
qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) + zqevap_ice(:,:,jl)* zmsk(:,:)
......@@ -2085,9 +2431,14 @@ CONTAINS
IF ( srcv(jpr_icb)%laction ) CALL iom_put('hflx_icb_cea' , - frcv(jpr_icb)%z3(:,:,1) * rLfus ) ! latent heat from icebergs melting
IF ( iom_use('hflx_rain_cea') ) & ! heat flux from rain (cell average)
& CALL iom_put('hflx_rain_cea' , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
IF ( iom_use('hflx_evap_cea') ) & ! heat flux from evap (cell average)
& CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) ) &
& * zcptn(:,:) * tmask(:,:,1) )
IF ( ln_couple_ocean_evap ) THEN
IF ( iom_use( 'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , frcv(jpr_tevp)%z3(:,:,1) * zcptn(:,:) * tmask(:,:,1) ) ! heat flux from evap (cell average)
ELSE
IF ( iom_use( 'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) &
& * zcptn(:,:) * tmask(:,:,1) ) ! heat flux from evap (cell average)
END IF
IF ( iom_use('hflx_prec_cea') ) & ! heat flux from all precip (cell avg)
& CALL iom_put('hflx_prec_cea' , sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
& + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
......@@ -2214,29 +2565,30 @@ CONTAINS
!
ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN !== conduction flux as surface forcing ==!
!
!! SELECT CASE( TRIM( sn_rcv_qtrice%cldes ) )
!! !
!! ! ! ===> here we receive the qtr_ice_top array from the coupler
!! CASE ('coupled')
!! IF (ln_scale_ice_flux) THEN
!! WHERE( a_i(:,:,:) > 1.e-10_wp )
!! zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
!! ELSEWHERE
!! zqtr_ice_top(:,:,:) = 0.0_wp
!! ENDWHERE
!! ELSE
!! zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:)
!! ENDIF
!!
!! ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
!! zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
!! zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
!!
!! ! if we are not getting this data from the coupler then assume zero (fully opaque ice)
!! CASE ('none')
zqtr_ice_top(:,:,:) = 0._wp
!! END SELECT
!
SELECT CASE( TRIM( sn_rcv_qtr%cldes ) )
!
! ! ===> here we receive the qtr_ice_top array from the coupler
CASE ('coupled')
IF (ln_scale_ice_flux) THEN
WHERE( a_i(:,:,:) > 0.0_wp ) zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:)) * a_i_last_couple(:,:,:) / a_i(:,:,:)
WHERE( a_i(:,:,:) <= 0.0_wp ) zqtr_ice_top(:,:,:) = 0.0_wp
ELSE
zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:))
ENDIF
! Add retrieved transmitted solar radiation onto the ice and total solar radiation
zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
! if we are not getting this data from the coupler then assume zero (fully opaque ice)
CASE ('none')
zqtr_ice_top(:,:,:) = 0._wp
CASE default
CALL ctl_stop( 'sbc_cpl_ice_flx: Invalid value for sn_rcv_qtr%cldes' )
END SELECT
!
ENDIF
IF( ln_mixcpl ) THEN
......@@ -2260,6 +2612,44 @@ CONTAINS
IF( ln_mixcpl ) THEN ; qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) + zqsr_oce(:,:)* zmsk(:,:)
ELSE ; qsr_oce(:,:) = zqsr_oce(:,:) ; ENDIF
IF( ln_mixcpl ) THEN
qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 ) ! total flux from blk
qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) + zqsr_tot(:,:)* zmsk(:,:)
DO jl = 1, jpl
qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) + zqsr_ice(:,:,jl)* zmsk(:,:)
END DO
ELSE
qsr_tot(:,: ) = zqsr_tot(:,: )
qsr_ice(:,:,:) = zqsr_ice(:,:,:)
ENDIF
! ! ========================= !
SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) ) ! d(qns)/dt !
! ! ========================= !
CASE ('coupled')
IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
ELSE
! Set all category values equal for the moment
DO jl=1,jpl
zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
ENDDO
ENDIF
CASE( 'none' )
zdqns_ice(:,:,:) = 0._wp
CASE default
CALL ctl_stop( 'sbc_cpl_ice_flx: Invalid value for sn_rcv_dqnsdt%cldes' )
END SELECT
IF( ln_mixcpl ) THEN
DO jl=1,jpl
dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
ENDDO
ELSE
dqns_ice(:,:,:) = zdqns_ice(:,:,:)
ENDIF
! ! ================== !
! ! ice skin temp. !
! ! ================== !
......@@ -2297,13 +2687,18 @@ CONTAINS
INTEGER, INTENT(in) :: kt
INTEGER, INTENT(in) :: Kbb, Kmm ! ocean model time level index
!
INTEGER :: ji, jj, jl ! dummy loop indices
INTEGER :: ji, jj, jl ! dummy loop indices
INTEGER :: ikchoix
INTEGER :: isec, info ! local integer
REAL(wp) :: zumax, zvmax
REAL(wp), DIMENSION(jpi,jpj) :: zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztmp3, ztmp4
REAL(wp), DIMENSION(jpi,jpj) :: ztmp5, ztmp6 ! RSRH temporary work arrays
! to avoid intent conflicts in repcmo calls
!!----------------------------------------------------------------------
!
IF (ln_timing) CALL timing_start('sbc_cpl_snd')
isec = ( kt - nit000 ) * NINT( rn_Dt ) ! date of exchanges
info = OASIS_idle
......@@ -2317,12 +2712,15 @@ CONTAINS
ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm) ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
ELSE
! we must send the surface potential temperature
IF( l_useCT ) THEN ; ztmp1(:,:) =eos_pt_from_ct( CASTSP(ts(:,:,1,jp_tem,Kmm)), CASTSP(ts(:,:,1,jp_sal,Kmm)) )
ELSE ; ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)
IF( l_useCT ) THEN
ztmp1(:,:) =eos_pt_from_ct( CASTSP(ts(:,:,1,jp_tem,Kmm)), CASTSP(ts(:,:,1,jp_sal,Kmm)) )
ELSE
ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)
ENDIF
!
SELECT CASE( sn_snd_temp%cldes)
CASE( 'oce only' ) ; ztmp1(:,:) = ztmp1(:,:) + rt0
CASE( 'oce only' )
ztmp1(:,:) = ztmp1(:,:) + rt0
CASE( 'oce and ice' ) ; ztmp1(:,:) = ztmp1(:,:) + rt0
SELECT CASE( sn_snd_temp%clcat )
CASE( 'yes' )
......@@ -2564,18 +2962,37 @@ CONTAINS
! ! ------------------------- !
! ! CO2 flux from PISCES !
! ! ------------------------- !
IF( ssnd(jps_co2)%laction .AND. l_co2cpl ) THEN
ztmp1(:,:) = oce_co2(:,:) * 1000. ! conversion in molC/m2/s
CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info )
IF( ssnd(jps_co2)%laction .AND. l_co2cpl ) THEN
ztmp1(:,:) = oce_co2(:,:) * 1000. ! conversion in molC/m2/s
CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info )
ENDIF
#if defined key_medusa
!
IF (ln_medusa) THEN
! ! ---------------------------------------------- !
! ! CO2 flux, DMS and chlorophyll from MEDUSA !
! ! ---------------------------------------------- !
IF ( ssnd(jps_bio_co2)%laction ) THEN
CALL cpl_snd( jps_bio_co2, isec, RESHAPE( CO2Flux_out_cpl, (/jpi,jpj,1/) ), info )
ENDIF
IF ( ssnd(jps_bio_dms)%laction ) THEN
CALL cpl_snd( jps_bio_dms, isec, RESHAPE( DMS_out_cpl, (/jpi,jpj,1/) ), info )
ENDIF
IF ( ssnd(jps_bio_chloro)%laction ) THEN
CALL cpl_snd( jps_bio_chloro, isec, RESHAPE( chloro_out_cpl, (/jpi,jpj,1/) ), info )
ENDIF
ENDIF
!
#endif
! ! ------------------------- !
IF( ssnd(jps_ocx1)%laction ) THEN ! Surface current !
! ! ------------------------- !
!
! j+1 j -----V---F
! surface velocity always sent from T point ! |
! j | T U
! [except for HadGEM3] j | T U
! | |
! j j-1 -I-------|
! (for I) | |
......@@ -2587,11 +3004,19 @@ CONTAINS
ELSE
SELECT CASE( TRIM( sn_snd_crt%cldes ) )
CASE( 'oce only' ) ! C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj ,1,Kmm) )
zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji ,jj-1,1,Kmm) )
END_2D
CASE( 'weighted oce and ice' ) ! Ocean and Ice on C-grid ==> T
IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj ,1,Kmm) )
zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji ,jj-1,1,Kmm) )
END_2D
ELSE
! Temporarily Changed for UKV
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = uu(ji,jj,1,Kmm)
zoty1(ji,jj) = vv(ji,jj,1,Kmm)
END_2D
ENDIF
CASE( 'weighted oce and ice' ) ! Ocean and Ice on C-grid ==> T
DO_2D( 0, 0, 0, 0 )
zotx1(ji,jj) = 0.5 * ( uu (ji,jj,1,Kmm) + uu (ji-1,jj ,1,Kmm) ) * zfr_l(ji,jj)
zoty1(ji,jj) = 0.5 * ( vv (ji,jj,1,Kmm) + vv (ji ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
......@@ -2607,6 +3032,7 @@ CONTAINS
& + 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj)
END_2D
END SELECT
CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
!
ENDIF
......@@ -2614,15 +3040,44 @@ CONTAINS
!
IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN ! Rotation of the components
! ! Ocean component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zotx1(:,:) = ztmp1(:,:) ! overwrite the components
zoty1(:,:) = ztmp2(:,:)
IF( ssnd(jps_ivx1)%laction ) THEN ! Ice component
CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zitx1(:,:) = ztmp1(:,:) ! overwrite the components
zity1(:,:) = ztmp2(:,:)
IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zotx1(:,:) = ztmp1(:,:) ! overwrite the components
zoty1(:,:) = ztmp2(:,:)
IF( ssnd(jps_ivx1)%laction ) THEN ! Ice component
CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component
CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component
zitx1(:,:) = ztmp1(:,:) ! overwrite the components
zity1(:,:) = ztmp2(:,:)
ENDIF
ELSE
! Temporary code for HadGEM3 - will be removed eventually.
! Only applies when we want uvel on U grid and vvel on V grid
! Rotate U and V onto geographic grid before sending.
DO_2D( 0, 0, 0, 0 )
ztmp1(ji,jj)=0.25*vmask(ji,jj,1) &
*(zotx1(ji,jj)+zotx1(ji-1,jj) &
+zotx1(ji,jj+1)+zotx1(ji-1,jj+1))
ztmp2(ji,jj)=0.25*umask(ji,jj,1) &
*(zoty1(ji,jj)+zoty1(ji+1,jj) &
+zoty1(ji,jj-1)+zoty1(ji+1,jj-1))
END_2D
ikchoix = -1
! zotx1 and zoty1 are input only to repcmo while ztmp5 and ztmp6
! are the newly calculated (output) values.
! Don't make the mistake of using zotx1 and zoty1 twice in this
! call for both input and output fields since it creates INTENT
! conflicts.
CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,ztmp5,ztmp6,ikchoix)
zotx1(:,:)=ztmp5(:,:)
zoty1(:,:)=ztmp6(:,:)
! Ensure any N fold and wrap columns are updated.
CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
ENDIF
ENDIF
!
......@@ -2783,6 +3238,7 @@ CONTAINS
IF( ssnd(jps_sstfrz)%laction ) CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info)
#endif
!
IF (ln_timing) CALL timing_stop('sbc_cpl_snd')
END SUBROUTINE sbc_cpl_snd
!!======================================================================
......
src/OCE/SBC/sbcrnf.F90
View file @ 4a9f4b18
......@@ -20,6 +20,7 @@ MODULE sbcrnf
USE sbc_oce ! surface boundary condition variables
USE eosbn2 ! Equation Of State
USE closea, ONLY: l_clo_rnf, clo_rnf ! closed seas
USE isf_oce
!
USE in_out_manager ! I/O manager
USE fldread ! read input field at current time step
......@@ -40,8 +41,10 @@ MODULE sbcrnf
LOGICAL :: ln_rnf_depth_ini !: depth river runoffs computed at the initialisation
REAL(wp) :: rn_rnf_max !: maximum value of the runoff climatologie (ln_rnf_depth_ini =T)
REAL(wp) :: rn_dep_max !: depth over which runoffs is spread (ln_rnf_depth_ini =T)
REAL(wp) :: tot_flux !: total iceberg flux (temporary variable)
INTEGER :: nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)
LOGICAL , PUBLIC :: ln_rnf_icb !: iceberg flux is specified in a file
LOGICAL :: ln_icb_mass !: Scale iceberg flux to match FW flux from coupled model
LOGICAL :: ln_rnf_tem !: temperature river runoffs attribute specified in a file
LOGICAL , PUBLIC :: ln_rnf_sal !: salinity river runoffs attribute specified in a file
TYPE(FLD_N) , PUBLIC :: sn_rnf !: information about the runoff file to be read
......@@ -119,8 +122,8 @@ CONTAINS
!
IF( .NOT. l_rnfcpl ) THEN
CALL fld_read ( kt, nn_fsbc, sf_rnf ) ! Read Runoffs data and provide it at kt ( runoffs + iceberg )
IF( ln_rnf_icb ) CALL fld_read ( kt, nn_fsbc, sf_i_rnf ) ! idem for iceberg flux if required
ENDIF
IF( ln_rnf_icb ) CALL fld_read ( kt, nn_fsbc, sf_i_rnf ) ! idem for iceberg flux if required
IF( ln_rnf_tem ) CALL fld_read ( kt, nn_fsbc, sf_t_rnf ) ! idem for runoffs temperature if required
IF( ln_rnf_sal ) CALL fld_read ( kt, nn_fsbc, sf_s_rnf ) ! idem for runoffs salinity if required
!
......@@ -138,14 +141,57 @@ CONTAINS
CALL iom_put( 'hflx_icb_cea' , -fwficb(:,:) * rLfus ) ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
ENDIF
ENDIF
!
IF( ln_rnf_icb ) THEN
fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1) ! updated runoff value at time step kt
IF( ln_icb_mass ) THEN
! Modify the Iceberg FW flux to be consistent with the change in
! mass of the Antarctic/Greenland ice sheet for FW conservation in
! coupled model. This isn't perfect as FW flux will go into ocean at
! wrong time of year but more important to maintain FW balance
tot_flux = SUM(fwficb(:,:)*e1e2t(:,:)*tmask_i(:,:)*greenland_icesheet_mask(:,:)) !Need to multiply by area to convert to kg/s
IF( lk_mpp ) CALL mpp_sum( 'icbclv', tot_flux )
IF( tot_flux > rsmall ) THEN
WHERE( greenland_icesheet_mask(:,:) == 1.0 ) &
& fwficb(:,:) = fwficb(:,:) * greenland_icesheet_mass_rate_of_change * rn_greenland_calving_fraction &
& / tot_flux
ELSE IF( rn_greenland_calving_fraction < rsmall ) THEN
WHERE( greenland_icesheet_mask(:,:) == 1.0 ) fwficb(:,:) = 0.0
ELSE
CALL CTL_STOP('STOP', 'No iceberg runoff data read in for Greenland. Check input file or set rn_greenland_calving_fraction=0.0')
ENDIF
tot_flux = SUM(fwficb(:,:)*e1e2t(:,:)*tmask_i(:,:)*antarctica_icesheet_mask(:,:))
IF( lk_mpp ) CALL mpp_sum( 'icbclv', tot_flux )
IF( tot_flux > rsmall ) THEN
WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
& fwficb(:,:) = fwficb(:,:) * antarctica_icesheet_mass_rate_of_change * rn_antarctica_calving_fraction &
& / (tot_flux + 1.0e-10_wp )
ELSE IF( rn_antarctica_calving_fraction < rsmall ) THEN
WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) fwficb(:,:) = 0.0
ELSE
CALL CTL_STOP('STOP', 'No iceberg runoff data read in for Antarctica. Check input file or set rn_antarctica_calving_fraction=0.0')
ENDIF
ENDIF
CALL lbc_lnk('sbcrnf',rnf,'T',1._wp,fwficb,'T',1._wp)
CALL iom_put( 'berg_melt' , fwficb(:,:) ) ! output iceberg flux
CALL iom_put( 'hflx_icb_cea' , fwficb(:,:) * rLfus ) ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
rnf(:,:) = rnf(:,:) + fwficb(:,:) ! fwficb isn't used anywhere else so add it to runoff here
qns_tot(:,:) = qns_tot(:,:) - fwficb(:,:) * rLfus ! TG: I think this is correct
ENDIF !
! ! set temperature & salinity content of runoffs
IF( ln_rnf_tem ) THEN ! use runoffs temperature data
rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0
CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp ) ! if missing data value use SST as runoffs temperature
rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rho0
rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rho0
END WHERE
WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )
! where fwf comes from melting of ice shelves or iceberg
rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rho0 - rnf(:,:) * rLfusisf * r1_rho0_rcp
END WHERE
ELSE ! use SST as runoffs temperature
!CEOD River is fresh water so must at least be 0 unless we consider ice
rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rho0
......@@ -158,11 +204,15 @@ CONTAINS
IF( iom_use('sflx_rnf_cea') ) CALL iom_put( 'sflx_rnf_cea', rnf_tsc(:,:,jp_sal) * rho0 ) ! output runoff salt flux (g/m2/s)
ENDIF
!
! Make sure rnf is up to date!
call lbc_lnk('rnf_div', rnf, 'T', 1.0)
! ! ---------------------------------------- !
IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 !
! ! ---------------------------------------- !
IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN !* Restart: read in restart file
IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields red in the restart file', lrxios
IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields read in the restart file', lrxios
CALL iom_get( numror, jpdom_auto, 'rnf_b' , rnf_b ) ! before runoff
CALL iom_get( numror, jpdom_auto, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem) ) ! before heat content of runoff
CALL iom_get( numror, jpdom_auto, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) ) ! before salinity content of runoff
......@@ -259,7 +309,7 @@ CONTAINS
REAL(wp) :: zacoef
REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl
!!
NAMELIST/namsbc_rnf/ cn_dir , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb, &
NAMELIST/namsbc_rnf/ cn_dir , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb, ln_icb_mass, &
& sn_rnf, sn_cnf , sn_i_rnf, sn_s_rnf , sn_t_rnf , sn_dep_rnf, &
& ln_rnf_mouth , rn_hrnf , rn_avt_rnf, rn_rfact, &
& ln_rnf_depth_ini , rn_dep_max , rn_rnf_max, nn_rnf_depth_file
......@@ -318,21 +368,22 @@ CONTAINS
IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) )
CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print )
!
IF( ln_rnf_icb ) THEN ! Create (if required) sf_i_rnf structure
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' iceberg flux read in a file'
ALLOCATE( sf_i_rnf(1), STAT=ierror )
IF( ierror > 0 ) THEN
CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' ) ; RETURN
ENDIF
ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1) )
IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) )
CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
ELSE
fwficb(:,:) = 0._wp
ENDIF
ENDIF
IF( ln_rnf_icb ) THEN ! Create (if required) sf_i_rnf structure
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' iceberg flux read in a file'
ALLOCATE( sf_i_rnf(1), STAT=ierror )
IF( ierror > 0 ) THEN
CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' ) ; RETURN
ENDIF
ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1) )
IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) )
CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
ELSE
fwficb(:,:) = 0._wp
ENDIF
!
IF( ln_rnf_tem ) THEN ! Create (if required) sf_t_rnf structure
IF(lwp) WRITE(numout,*)
......
src/OCE/SBC/sbcssm.F90
View file @ 4a9f4b18
......@@ -59,13 +59,28 @@ CONTAINS
INTEGER :: ji, jj ! loop index
REAL(wp) :: zcoef, zf_sbc ! local scalar
REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts
!!---------------------------------------------------------------------
CHARACTER(len=4),SAVE :: stype
!!---------------------------------------------------------------------
IF( kt == nit000 ) THEN
IF( ln_TEOS10 ) THEN
stype='abs' ! teos-10: using absolute salinity (sst is converted to potential temperature for the surface module)
ELSE IF( ln_SEOS ) THEN
stype='seos' ! seos: using Simplified Equation of state (sst is converted to potential temperature for the surface module)
ELSE
stype='pra' ! eos-80: using practical salinity
ENDIF
ENDIF
!
! !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
zts(:,:,jp_tem) = ts(:,:,1,jp_tem,Kmm)
zts(:,:,jp_sal) = ts(:,:,1,jp_sal,Kmm)
!
! ! ---------------------------------------- !
! !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain
! ! otherwise restartability and reproducibility are broken
! ! computed in traqsr only on the inner domain
CALL lbc_lnk( 'sbc_ssm', fraqsr_1lev(:,:), 'T', 1._wp )
!
! ! ---------------------------------------- !
IF( nn_fsbc == 1 ) THEN ! Instantaneous surface fields !
! ! ---------------------------------------- !
ssu_m(:,:) = uu(:,:,1,Kbb)
......@@ -171,8 +186,8 @@ CONTAINS
IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN ! Mean value at each nn_fsbc time-step !
CALL iom_put( 'ssu_m', ssu_m )
CALL iom_put( 'ssv_m', ssv_m )
CALL iom_put( 'sst_m', sst_m )
CALL iom_put( 'sss_m', sss_m )
CALL iom_put( 'sst_m_pot', sst_m )
CALL iom_put( 'sss_m_'//stype, sss_m )
CALL iom_put( 'ssh_m', ssh_m )
CALL iom_put( 'e3t_m', e3t_m )
CALL iom_put( 'frq_m', frq_m )
......
src/OCE/TRA/eosbn2.F90
View file @ 4a9f4b18
......@@ -45,6 +45,7 @@ MODULE eosbn2
!
USE in_out_manager ! I/O manager
USE lib_mpp ! for ctl_stop
USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
USE prtctl ! Print control
USE timing ! Timing
......
src/OCE/TRA/traadv_fct.F90
View file @ 4a9f4b18
......@@ -87,6 +87,7 @@ CONTAINS
REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - -
REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - -
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwz, ztu, ztv, zltu, zltv, ztw
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi_in, ztw_in !temp arrays to avoid intent conflicts
REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zwx, zwy
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup
......@@ -194,7 +195,9 @@ CONTAINS
END_3D
IF ( ll_zAimp ) THEN
CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
! We need to use separate copies of zwi to avoid intent conflicts!
zwi_in(:,:,:) = zwi(:,:,:)
CALL tridia_solver( zwdia, zwsup, zwinf, zwi_in, zwi , 0 )
!
ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)
......@@ -314,7 +317,10 @@ CONTAINS
ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
END_3D
!
CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
! We need to use separate copies of ztw to avoid intent conflicts!
ztw_in(:,:,:) = ztw(:,:,:)
CALL tridia_solver( zwdia, zwsup, zwinf, ztw_in, ztw , 0 )
!
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)
zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
......@@ -724,6 +730,8 @@ CONTAINS
REAL(wp) :: zwy_jm1, zfm_ui, zfm_ui_m1, zfm_vj, zfm_vj_m1, zfm_wk, zC2t_v, zC4t_v ! - -
REAL(wp) :: ztu, ztv, ztu_im1, ztu_ip1, ztv_jm1, ztv_jp1
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi_in, ztw_in ! Read-only copies to avoid INTENT conflicts
! in calls to tridia_solver
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup
LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection
......@@ -739,6 +747,8 @@ CONTAINS
zwy_3d(:,:,:) = 0._wp
zwz(:,:,:) = 0._wp
zwi(:,:,:) = 0._wp
zwt(:,:,:) = 0._wp
!
l_trd = .FALSE. ! set local switches
l_hst = .FALSE.
......@@ -820,7 +830,10 @@ CONTAINS
END DO
IF ( ll_zAimp ) THEN
CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 )
! We need to use separate copies of zwi to avoid intent conflicts!
zwi_in(:,:,:) = zwi(:,:,:)
CALL tridia_solver( zwdia, zwsup, zwinf, zwi_in, zwi , 0 )
!
ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)
......@@ -934,7 +947,9 @@ CONTAINS
ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
END_3D
!
CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
! We need to use separate copies of ztw to avoid intent conflicts!
ztw_in(:,:,:) = ztw(:,:,:)
CALL tridia_solver( zwdia, zwsup, zwinf, ztw_in, ztw , 0 )
!
DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)
zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
......
src/OCE/TRA/traqsr.F90
View file @ 4a9f4b18
......@@ -48,6 +48,7 @@ MODULE traqsr
LOGICAL , PUBLIC :: ln_qsr_2bd !: 2 band light absorption flag
LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag
INTEGER , PUBLIC :: nn_chldta !: use Chlorophyll data (=1) or not (=0)
REAL(wp), PUBLIC :: rn_chl_conc !: Chlorophyll concentration (for nn_chldta=0)
REAL(wp), PUBLIC :: rn_abs !: fraction absorbed in the very near surface (RGB & 2 bands)
REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands)
REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands)
......@@ -221,7 +222,7 @@ CONTAINS
ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
END_3D
ELSE !* constant chlorophyll
zchl = 0.05
zchl = rn_chl_conc
! NB. make sure constant value is such that:
zchl = MIN( 10. , MAX( 0.03, zchl ) )
! Convert chlorophyll value to attenuation coefficient look-up table index
......@@ -341,14 +342,14 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ios, irgb, ierror, ioptio ! local integer
REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars
REAL(wp) :: zz1, zc2 , zc3, zchl ! - -
REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars
REAL(wp) :: zz1, zc2 , zc3, zchl ! - -
!
CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read
!!
NAMELIST/namtra_qsr/ sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, &
& nn_chldta, rn_abs, rn_si0, rn_si1
& nn_chldta, rn_chl_conc, rn_abs, rn_si0, rn_si1
!!----------------------------------------------------------------------
!
READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
......@@ -367,6 +368,7 @@ CONTAINS
WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd
WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio
WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta
WRITE(numout,*) ' Chlorophyll concentration (for nn_chldta=0) rn_chl_conc = ', rn_chl_conc
WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs
WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0
WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1
......@@ -415,7 +417,7 @@ CONTAINS
& 'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
ENDIF
IF( nqsr == np_RGB ) THEN ! constant Chl
IF(lwp) WRITE(numout,*) ' ==>>> Constant Chlorophyll concentration = 0.05'
IF(lwp) WRITE(numout,*) ' ==>>> Constant Chlorophyll concentration = ', rn_chl_conc
ENDIF
!
CASE( np_2BD ) !== 2 bands light penetration ==!
......
src/OCE/ZDF/zdf_oce.F90
View file @ 4a9f4b18
......@@ -39,6 +39,7 @@ MODULE zdf_oce
! ! gravity wave-induced vertical mixing
LOGICAL , PUBLIC :: ln_zdfswm !: surface wave-induced mixing flag
LOGICAL , PUBLIC :: ln_zdfiwm !: internal wave-induced mixing flag
LOGICAL , PUBLIC :: ln_zdftmx !: old tidal mixing scheme (Simmons et al 2004)
LOGICAL , PUBLIC :: ln_zdfmfc !: convection: eddy diffusivity Mass Flux Convection
! ! coefficients
REAL(wp), PUBLIC :: rn_avm0 !: vertical eddy viscosity (m2/s)
......
src/OCE/ZDF/zdfmxl.F90
View file @ 4a9f4b18
......@@ -15,6 +15,7 @@ MODULE zdfmxl
USE dom_oce ! ocean space and time domain variables
USE trc_oce , ONLY: l_offline ! ocean space and time domain variables
USE zdf_oce ! ocean vertical physics
USE eosbn2 ! for zdf_mxl_zint
!
USE in_out_manager ! I/O manager
USE prtctl ! Print control
......@@ -26,15 +27,27 @@ MODULE zdfmxl
PRIVATE
PUBLIC zdf_mxl, zdf_mxl_turb, zdf_mxl_alloc ! called by zdfphy.F90
PUBLIC zdf_mxl_zint ! called by diahth.F90
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nmln !: number of level in the mixed layer (used by LDF, ZDF, TRD, TOP)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld !: mixing layer depth (turbocline) [m] (used by TOP)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] (used by LDF)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: depth of the last T-point inside the mixed layer [m] (used by LDF)
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hmld_zint !: vertically-interpolated mixed layer depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: htc_mld ! Heat content of hmld_zint
LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ll_found ! Is T_b to be found by interpolation ?
LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ll_belowml ! Flag points below mixed layer when ll_found=F
REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth
REAL(wp), PUBLIC :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth
TYPE, PUBLIC :: MXL_ZINT !: Structure for MLD defs
INTEGER :: mld_type ! mixed layer type
REAL(wp) :: zref ! depth of initial T_ref
REAL(wp) :: dT_crit ! Critical temp diff
REAL(wp) :: iso_frac ! Fraction of rn_dT_crit
END TYPE MXL_ZINT
!! * Substitutions
# include "do_loop_substitute.h90"
# include "single_precision_substitute.h90"
......@@ -52,7 +65,8 @@ CONTAINS
!!----------------------------------------------------------------------
zdf_mxl_alloc = 0 ! set to zero if no array to be allocated
IF( .NOT. ALLOCATED( nmln ) ) THEN
ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc )
ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj), &
& htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
!
CALL mpp_sum ( 'zdfmxl', zdf_mxl_alloc )
IF( zdf_mxl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_alloc: failed to allocate arrays.' )
......@@ -86,6 +100,8 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
IF(lwp) WRITE(numout,*) '~~~~~~~ '
! ! allocate zdfmxl arrays
IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
ENDIF
ENDIF
!
......@@ -120,6 +136,352 @@ CONTAINS
!
END SUBROUTINE zdf_mxl
SUBROUTINE zdf_mxl_zint_mld( sf , Kmm)
!!----------------------------------------------------------------------------------
!! *** ROUTINE zdf_mxl_zint_mld ***
!
! Calculate vertically-interpolated mixed layer depth diagnostic.
!
! This routine can calculate the mixed layer depth diagnostic suggested by
! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate
! vertically-interpolated mixed-layer depth diagnostics with other parameter
! settings set in the namzdf_mldzint namelist.
!
! If mld_type=1 the mixed layer depth is calculated as the depth at which the
! density has increased by an amount equivalent to a temperature difference of
! 0.8C at the surface.
!
! For other values of mld_type the mixed layer is calculated as the depth at
! which the temperature differs by 0.8C from the surface temperature.
!
! David Acreman, Daley Calvert
!
!!-----------------------------------------------------------------------------------
TYPE(MXL_ZINT), INTENT(in) :: sf
INTEGER, INTENT(in) :: Kmm ! ocean time level index
! Diagnostic criteria
INTEGER :: nn_mld_type ! mixed layer type
REAL(wp) :: rn_zref ! depth of initial T_ref
REAL(wp) :: rn_dT_crit ! Critical temp diff
REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used
! Local variables
REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value
INTEGER, DIMENSION(jpi,jpj) :: ikmt ! number of active tracer levels
INTEGER, DIMENSION(jpi,jpj) :: ik_ref ! index of reference level
INTEGER, DIMENSION(jpi,jpj) :: ik_iso ! index of last uniform temp level
REAL, DIMENSION(jpi,jpj,jpk) :: zT ! Temperature or density
REAL, DIMENSION(jpi,jpj) :: ppzdep ! depth for use in calculating d(rho)
REAL, DIMENSION(jpi,jpj) :: zT_ref ! reference temperature
REAL :: zT_b ! base temperature
REAL, DIMENSION(jpi,jpj,jpk) :: zdTdz ! gradient of zT
REAL, DIMENSION(jpi,jpj,jpk) :: zmoddT ! Absolute temperature difference
REAL :: zdz ! depth difference
REAL :: zdT ! temperature difference
REAL, DIMENSION(jpi,jpj) :: zdelta_T ! difference critereon
REAL, DIMENSION(jpi,jpj) :: zRHO1, zRHO2 ! Densities
INTEGER :: ji, jj, jk ! loop counter
!!-------------------------------------------------------------------------------------
!
! Unpack structure
nn_mld_type = sf%mld_type
rn_zref = sf%zref
rn_dT_crit = sf%dT_crit
rn_iso_frac = sf%iso_frac
! Set the mixed layer depth criterion at each grid point
IF( nn_mld_type == 0 ) THEN
zdelta_T(:,:) = rn_dT_crit
zT(:,:,:) = rhop(:,:,:)
ELSE IF( nn_mld_type == 1 ) THEN
ppzdep(:,:)=0.0
call eos ( ts(:,:,1,:,Kmm), ppzdep(:,:), zRHO1(:,:) )
! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST
! [assumes number of tracers less than number of vertical levels]
zT(:,:,1:jpts)=ts(:,:,1,1:jpts,Kmm)
zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit
CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )
zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rho0
! RHO from eos (2d version) doesn't calculate north or east halo:
CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. )
zT(:,:,:) = rhop(:,:,:)
ELSE
zdelta_T(:,:) = rn_dT_crit
zT(:,:,:) = ts(:,:,:,jp_tem,Kmm)
END IF
! Calculate the gradient of zT and absolute difference for use later
DO jk = 1 ,jpk-2
zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w(:,:,jk+1,Kmm)
zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )
END DO
! Find density/temperature at the reference level (Kara et al use 10m).
! ik_ref is the index of the box centre immediately above or at the reference level
! Find rn_zref in the array of model level depths and find the ref
! density/temperature by linear interpolation.
DO jk = jpkm1, 2, -1
WHERE ( gdept(:,:,jk,Kmm) > rn_zref )
ik_ref(:,:) = jk - 1
zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept(:,:,jk-1,Kmm) )
END WHERE
END DO
! If the first grid box centre is below the reference level then use the
! top model level to get zT_ref
WHERE ( gdept(:,:,1,Kmm) > rn_zref )
zT_ref = zT(:,:,1)
ik_ref = 1
END WHERE
! The number of active tracer levels is 1 less than the number of active w levels
ikmt(:,:) = mbkt(:,:) - 1
! Initialize / reset
ll_found(:,:) = .false.
IF ( rn_iso_frac - zepsilon > 0. ) THEN
! Search for a uniform density/temperature region where adjacent levels
! differ by less than rn_iso_frac * deltaT.
! ik_iso is the index of the last level in the uniform layer
! ll_found indicates whether the mixed layer depth can be found by interpolation
ik_iso(:,:) = ik_ref(:,:)
DO jj = 1, jpj ! Changed from nlcj
DO ji = 1, jpi ! Changed from nlci
!CDIR NOVECTOR
DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1
IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN
ik_iso(ji,jj) = jk
ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )
EXIT
END IF
END DO
END DO
END DO
! Use linear interpolation to find depth of mixed layer base where possible
hmld_zint(:,:) = rn_zref
DO jj = 1, jpj
DO ji = 1, jpi
IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN
zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )
hmld_zint(ji,jj) = gdept(ji,jj,ik_iso(ji,jj),Kmm) + zdz
END IF
END DO
END DO
END IF
! If ll_found = .false. then calculate MLD using difference of zdelta_T
! from the reference density/temperature
! Prevent this section from working on land points
WHERE ( tmask(:,:,1) /= 1.0 )
ll_found = .true.
END WHERE
DO jk=1, jpk
ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)
END DO
! Set default value where interpolation cannot be used (ll_found=false)
DO jj = 1, jpj
DO ji = 1, jpi
IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = gdept(ji,jj,ikmt(ji,jj),Kmm)
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
!CDIR NOVECTOR
DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)
IF ( ll_found(ji,jj) ) EXIT
IF ( ll_belowml(ji,jj,jk) ) THEN
zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )
zdT = zT_b - zT(ji,jj,jk-1)
zdz = zdT / zdTdz(ji,jj,jk-1)
hmld_zint(ji,jj) = gdept(ji,jj,jk-1,Kmm) + zdz
EXIT
END IF
END DO
END DO
END DO
hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)
!
END SUBROUTINE zdf_mxl_zint_mld
SUBROUTINE zdf_mxl_zint_htc( kt , Kmm)
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_mxl_zint_htc ***
!!
!! ** Purpose :
!!
!! ** Method :
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step index
INTEGER, INTENT(in) :: Kmm ! ocean time level index
INTEGER :: ji, jj, jk
INTEGER :: ikmax
REAL(wp) :: zc, zcoef
!
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilevel
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zthick_0, zthick
!!----------------------------------------------------------------------
IF( .NOT. ALLOCATED(ilevel) ) THEN
ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), &
& zthick(jpi,jpj), STAT=ji )
IF( lk_mpp ) CALL mpp_sum( 'zdfmxl', ji )
IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' )
ENDIF
! Find last whole model T level above the MLD
ilevel(:,:) = 0
zthick_0(:,:) = 0._wp
DO jk = 1, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
zthick_0(ji,jj) = zthick_0(ji,jj) + e3t(ji,jj,jk,Kmm)
IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk
END DO
END DO
WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2)
WRITE(numout,*) 'gdepw(jk+1 =',jk+1,') =',gdepw(2,2,jk+1,Kmm)
END DO
! Surface boundary condition
IF( ln_linssh ) THEN ; zthick(:,:) = ssh(:,:,Kmm) ; htc_mld(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) * tmask(:,:,1)
ELSE ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp
ENDIF
! Deepest whole T level above the MLD
ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 )
! Integration down to last whole model T level
DO jk = 1, ikmax
DO jj = 1, jpj
DO ji = 1, jpi
zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1 ) ) ! 0 below ilevel
zthick(ji,jj) = zthick(ji,jj) + zc
htc_mld(ji,jj) = htc_mld(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
END DO
END DO
END DO
! Subsequent partial T level
zthick(:,:) = hmld_zint(:,:) - zthick(:,:) ! remaining thickness to reach MLD
DO jj = 1, jpj
DO ji = 1, jpi
htc_mld(ji,jj) = htc_mld(ji,jj) + ts(ji,jj,ilevel(ji,jj)+1,jp_tem,Kmm) &
& * MIN( e3t(ji,jj,ilevel(ji,jj)+1,Kmm), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)
END DO
END DO
WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2)
! Convert to heat content
zcoef = rho0 * rcp
htc_mld(:,:) = zcoef * htc_mld(:,:)
END SUBROUTINE zdf_mxl_zint_htc
SUBROUTINE zdf_mxl_zint( kt , Kmm)
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_mxl_zint ***
!!
!! ** Purpose :
!!
!! ** Method :
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time-step index
INTEGER, INTENT(in) :: Kmm ! ocean time level index
INTEGER :: ios
INTEGER :: jn
INTEGER :: nn_mld_diag = 0 ! number of diagnostics
CHARACTER(len=1) :: cmld
LOGICAL, SAVE, DIMENSION(5) :: l_mld, l_htc
TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
TYPE(MXL_ZINT), SAVE, DIMENSION(5) :: mld_diags
NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
!!----------------------------------------------------------------------
IF( kt == nit000 ) THEN
READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' )
READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' )
IF(lwm) WRITE ( numond, namzdf_mldzint )
IF( nn_mld_diag > 5 ) CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' )
l_mld(:) = .FALSE. ; l_htc(:) = .FALSE.
mld_diags(1) = sn_mld1
mld_diags(2) = sn_mld2
mld_diags(3) = sn_mld3
mld_diags(4) = sn_mld4
mld_diags(5) = sn_mld5
IF( nn_mld_diag > 0 ) THEN
IF( lwp ) THEN
WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================'
WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)'
ENDIF
DO jn = 1, nn_mld_diag
! Whether the diagnostic is requested
WRITE(cmld,'(I1)') jn
IF( iom_use( "mldzint_"//cmld ) ) l_mld(jn) = .TRUE.
IF( iom_use( "mldhtc_"//cmld ) ) l_htc(jn) = .TRUE.
IF( lwp ) THEN
WRITE(numout,*) 'MLD criterion',jn,':'
WRITE(numout,*) ' nn_mld_type =', mld_diags(jn)%mld_type
WRITE(numout,*) ' rn_zref =' , mld_diags(jn)%zref
WRITE(numout,*) ' rn_dT_crit =' , mld_diags(jn)%dT_crit
WRITE(numout,*) ' rn_iso_frac =', mld_diags(jn)%iso_frac
ENDIF
END DO
WRITE(numout,*) '===================================================================='
ENDIF
ENDIF
IF( nn_mld_diag > 0 ) THEN
DO jn = 1, nn_mld_diag
WRITE(cmld,'(I1)') jn
IF( l_mld(jn) .OR. l_htc(jn) ) THEN
CALL zdf_mxl_zint_mld( mld_diags(jn), Kmm)
IF( l_mld(jn) ) CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) )
IF( l_htc(jn) ) THEN
CALL zdf_mxl_zint_htc( kt, Kmm )
CALL iom_put( "mldhtc_"//cmld, htc_mld(:,:) )
ENDIF
ENDIF
END DO
ENDIF
END SUBROUTINE zdf_mxl_zint
SUBROUTINE zdf_mxl_turb( kt, Kmm )
!!----------------------------------------------------------------------
......
src/OCE/ZDF/zdfphy.F90
View file @ 4a9f4b18
......@@ -24,6 +24,7 @@ MODULE zdfphy
USE zdfevd ! vertical physics: convection via enhanced vertical diffusion
USE zdfmfc ! vertical physics: Mass Flux Convection
USE zdfiwm ! vertical physics: internal wave-induced mixing
USE zdftmx ! vertical physics: old tidal mixing scheme (Simmons et al 2004)
USE zdfswm ! vertical physics: surface wave-induced mixing
USE zdfmxl ! vertical physics: mixed layer
USE tranpc ! convection: non penetrative adjustment
......@@ -89,6 +90,7 @@ CONTAINS
& ln_zdfnpc, nn_npc , nn_npcp, & ! convection : npc
& ln_zdfddm, rn_avts, rn_hsbfr, & ! double diffusion
& ln_zdfswm, & ! surface wave-induced mixing
& ln_zdftmx, & ! old tidal mixing scheme (Simmons et al 2004)
& ln_zdfiwm, & ! internal - - -
& ln_zad_Aimp, & ! apdative-implicit vertical advection
& rn_avm0, rn_avt0, nn_avb, nn_havtb ! coefficients
......@@ -227,6 +229,7 @@ CONTAINS
!
! !== gravity wave-driven mixing ==!
IF( ln_zdfiwm ) CALL zdf_iwm_init ! internal wave-driven mixing
IF( ln_zdftmx ) CALL zdf_tmx_init ! old tidal mixing scheme (Simmons et al)
IF( ln_zdfswm ) CALL zdf_swm_init ! surface wave-driven mixing
! !== top/bottom friction ==!
......@@ -304,7 +307,7 @@ CONTAINS
ENDIF
!
SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points
CASE( np_RIC ) ; CALL zdf_ric( kt, Kmm, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz
CASE( np_RIC ) ; CALL zdf_ric( kt, Kbb, Kmm, avm_k, avt_k ) ! Richardson number dependent Kz
CASE( np_TKE ) ; CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz
CASE( np_GLS ) ; CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz
CASE( np_OSM ) ; CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k ) ! OSMOSIS closure scheme for Kz
......@@ -348,6 +351,7 @@ CONTAINS
! !* wave-induced mixing
IF( ln_zdfswm ) CALL zdf_swm( kt, Kmm, avm, avt, avs ) ! surface wave (Qiao et al. 2004)
IF( ln_zdfiwm ) CALL zdf_iwm( kt, Kmm, avm, avt, avs ) ! internal wave (de Lavergne et al 2017)
IF( ln_zdftmx ) CALL zdf_tmx( kt, Kmm, avm, avt, avs ) ! old tidal mixing scheme (Simmons et al 2004)
! !* Lateral boundary conditions (sign unchanged)
IF(nn_hls==1) THEN
......
src/OCE/ZDF/zdfric.F90
View file @ 4a9f4b18
......@@ -106,7 +106,7 @@ CONTAINS
END SUBROUTINE zdf_ric_init
SUBROUTINE zdf_ric( kt, Kmm, p_sh2, p_avm, p_avt )
SUBROUTINE zdf_ric( kt, Kbb, Kmm, p_avm, p_avt )
!!----------------------------------------------------------------------
!! *** ROUTINE zdfric ***
!!
......@@ -146,18 +146,22 @@ CONTAINS
!! PFJ Lermusiaux 2001.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term
INTEGER , INTENT(in ) :: Kmm , Kbb ! ocean time level index
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zcfRi, zav, zustar, zhek ! local scalars
REAL(wp) :: zcfRi, zav, zustar, zhek, zdku, zdkv, zwx ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zh_ekm ! 2D workspace
!!----------------------------------------------------------------------
!
! !== avm and avt = F(Richardson number) ==!
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! coefficient = F(richardson number) (avm-weighted Ri)
zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) ) )
zdku = 0.5 / e3uw(ji,jj,jk,Kbb) * ( uu(ji-1,jj,jk-1,Kbb) + uu(ji,jj,jk-1,Kbb) &
& - uu(ji-1,jj,jk,Kbb) - uu(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk)
zdkv = 0.5 / e3vw(ji,jj,jk,Kbb) * ( vv(ji,jj-1,jk-1,Kbb) + vv(ji,jj,jk-1,Kbb) &
& - vv(ji,jj-1,jk,Kbb) - vv(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk)
zwx = zdku*zdku + zdkv*zdkv
zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , rn2(ji,jj,jk) / ( zwx + 1.e-20 ) ) )
zav = rn_avmri * zcfRi**nn_ric
! ! avm and avt coefficients
p_avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk)
......
src/OCE/ZDF/zdftke.F90
View file @ 4a9f4b18
......@@ -84,6 +84,7 @@ MODULE zdftke
REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3)
INTEGER :: nn_htau ! type of tke profile of penetration (=0/1)
REAL(wp) :: rn_htau_scaling ! a acaling factor to apply to the penetration of TKE
INTEGER :: nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
INTEGER :: nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean
......@@ -96,7 +97,7 @@ MODULE zdftke
REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3)
REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation
......@@ -296,17 +297,13 @@ CONTAINS
! ! 1/2 (W_lc)^2 = MAX( u* u_s + v* v_s , 0 ) only the positive part
!!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s
!!gm ! so we will overestimate the LC velocity.... !!gm I will do the work if !LC have an effect !
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
!!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) )
zWlc2(ji,jj) = 0.5_wp * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )
END_2D
!
! Projection of Stokes drift in the wind stress direction
!
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
ztauj = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1)
z1_norm = 1._wp / MAX( ztaui*ztaui+ztauj*ztauj, 1.e-12 ) * tmask(ji,jj,1)
zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2
END_2D
ELSE ! Surface Stokes drift deduced from surface stress
......@@ -737,6 +734,7 @@ CONTAINS
& rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &
& rn_mxl0 , nn_mxlice, rn_mxlice, &
& nn_pdl , ln_lc , rn_lc , &
& rn_htau_scaling , &
& nn_etau , nn_htau , rn_efr , nn_eice , &
& nn_bc_surf, nn_bc_bot, ln_mxhsw
!!----------------------------------------------------------------------
......@@ -793,6 +791,7 @@ CONTAINS
ENDIF
WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau
WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau
WRITE(numout,*) ' scaling factor for tke penetration depth rn_htau_scaling = ', rn_htau_scaling
WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr
WRITE(numout,*) ' langmuir & surface wave breaking under ice nn_eice = ', nn_eice
SELECT CASE( nn_eice )
......@@ -823,7 +822,7 @@ CONTAINS
! !* Check of some namelist values
IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1, 2 or 3' )
IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1' )
IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' )
IF( ( nn_htau < 0 .OR. nn_htau > 1 ) .AND. nn_htau .NE. 4 .AND. nn_htau .NE. 5 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 , 4 or 5 ' )
IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
!
IF( ln_mxl0 ) THEN
......@@ -835,9 +834,24 @@ CONTAINS
IF( nn_etau /= 0 ) THEN
SELECT CASE( nn_htau ) ! Choice of the depth of penetration
CASE( 0 ) ! constant depth penetration (here 10 meters)
htau(:,:) = 10._wp
htau(:,:) = rn_htau_scaling*10._wp
CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees
htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) )
CASE( 4 ) ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
IF( gphit(ji,jj) <= 0._wp ) THEN
htau(ji,jj) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* rn_htau_scaling*ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) )
ELSE
htau(ji,jj) = MAX( 0.5_wp, MIN( 10._wp, 45._wp* rn_htau_scaling*ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) )
ENDIF
END_2D
CASE( 5 ) ! Variation on case 4 with a steeper ramp further south in Southern Hemisphere
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
htau(ji,jj) = MAX( 0.5_wp, MIN( 10._wp, 45._wp* rn_htau_scaling*ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) )
IF( gphit(ji,jj) <= -40._wp ) THEN
htau(ji,jj) = htau(ji,jj) + MIN( 20._wp, 135._wp * ABS( SIN( rpi/180._wp * (gphit(ji,jj) + 40.0) ) ) )
ENDIF
END_2D
END SELECT
ENDIF
! !* read or initialize all required files
......
src/OCE/ZDF/zdftmx.F90 0 → 100644
View file @ 4a9f4b18
MODULE zdftmx
!!========================================================================
!! *** MODULE zdftmx ***
!! Ocean physics: vertical tidal mixing coefficient
!!========================================================================
!! History : 1.0 ! 2004-04 (L. Bessieres, G. Madec) Original code
!! - ! 2006-08 (A. Koch-Larrouy) Indonesian strait
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! 'key_zdftmx' Tidal vertical mixing
!!----------------------------------------------------------------------
!! zdf_tmx : global momentum & tracer Kz with tidal induced Kz
!! tmx_itf : Indonesian momentum & tracer Kz with tidal induced Kz
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE dom_oce ! ocean space and time domain variables
USE zdf_oce ! ocean vertical physics variables
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE eosbn2 ! ocean equation of state
USE phycst ! physical constants
USE prtctl ! Print control
USE in_out_manager ! I/O manager
USE iom ! I/O Manager
USE lib_mpp ! MPP library
USE timing ! Timing
USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
IMPLICIT NONE
PRIVATE
PUBLIC zdf_tmx ! called in step module
PUBLIC zdf_tmx_init ! called in opa module
PUBLIC zdf_tmx_alloc ! called in nemogcm module
LOGICAL, PUBLIC, PARAMETER :: lk_zdftmx = .TRUE. !: tidal mixing flag
! !!* Namelist namzdf_tmx : tidal mixing *
REAL(wp) :: rn_htmx ! vertical decay scale for turbulence (meters)
REAL(wp) :: rn_n2min ! threshold of the Brunt-Vaisala frequency (s-1)
REAL(wp) :: rn_tfe ! tidal dissipation efficiency (St Laurent et al. 2002)
REAL(wp) :: rn_me ! mixing efficiency (Osborn 1980)
LOGICAL :: ln_tmx_itf ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
REAL(wp) :: rn_tfe_itf ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: en_tmx ! energy available for tidal mixing (W/m2)
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mask_itf ! mask to use over Indonesian area
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: az_tmx ! coefficient used to evaluate the tidal induced Kz
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OPA 4.0 , NEMO Consortium (2011)
!! $Id: zdftmx.F90 8788 2017-11-22 18:01:02Z davestorkey $
!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
INTEGER FUNCTION zdf_tmx_alloc()
!!----------------------------------------------------------------------
!! *** FUNCTION zdf_tmx_alloc ***
!!----------------------------------------------------------------------
ALLOCATE(en_tmx(jpi,jpj), mask_itf(jpi,jpj), az_tmx(jpi,jpj,jpk), STAT=zdf_tmx_alloc )
!
IF( lk_mpp ) CALL mpp_sum ('zdftmx', zdf_tmx_alloc )
IF( zdf_tmx_alloc /= 0 ) CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
END FUNCTION zdf_tmx_alloc
SUBROUTINE zdf_tmx( kt, Kmm, p_avm, p_avt, p_avs)
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_tmx ***
!!
!! ** Purpose : add to the vertical mixing coefficients the effect of
!! tidal mixing (Simmons et al 2004).
!!
!! ** Method : - tidal-induced vertical mixing is given by:
!! Kz_tides = az_tmx / max( rn_n2min, N^2 )
!! where az_tmx is a coefficient that specified the 3D space
!! distribution of the faction of tidal energy taht is used
!! for mixing. Its expression is set in zdf_tmx_init routine,
!! following Simmons et al. 2004.
!! NB: a specific bounding procedure is performed on av_tide
!! so that the input tidal energy is actually almost used. The
!! basic maximum value is 60 cm2/s, but values of 300 cm2/s
!! can be reached in area where bottom stratification is too
!! weak.
!!
!! - update av_tide in the Indonesian Through Flow area
!! following Koch-Larrouy et al. (2007) parameterisation
!! (see tmx_itf routine).
!!
!! - update the model vertical eddy viscosity and diffusivity:
!! avt = avt + av_tides
!! avm = avm + av_tides
!!
!! ** Action : avt, avm increased by tidal mixing
!!
!! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
!! Koch-Larrouy et al. 2007, GRL.
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt , Kmm ! ocean time-step , time level
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm ! momentum Kz (w-points)
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avt, p_avs ! tracer Kz (w-points)
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: ztpc ! scalar workspace
REAL(wp), DIMENSION(jpi,jpj) :: zkz
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zav_tide
!!----------------------------------------------------------------------
!
! ! ----------------------- !
! ! Standard tidal mixing ! (compute zav_tide)
! ! ----------------------- !
! !* First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
zav_tide(:,:,:) = MIN( 60.e-4, az_tmx(:,:,:) / MAX( rn_n2min, rn2(:,:,:) ) )
zkz(:,:) = 0.e0 !* Associated potential energy consummed over the whole water column
DO jk = 2, jpkm1
zkz(:,:) = zkz(:,:) + e3w(:,:,jk,Kmm) * MAX( 0.e0, rn2(:,:,jk) ) * rho0 * zav_tide(:,:,jk) * wmask(:,:,jk)
END DO
DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
DO ji = 1, jpi
IF( zkz(ji,jj) /= 0.e0 ) zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
END DO
END DO
DO jk = 2, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
DO ji = 1, jpi
zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s
END DO
END DO
END DO
IF( kt == nit000 ) THEN !* check at first time-step: diagnose the energy consumed by zav_tide
ztpc = 0.e0
DO jk= 1, jpk
DO jj= 1, jpj
DO ji= 1, jpi
ztpc = ztpc + e3w(ji,jj,jk,Kmm) * e1t(ji,jj) * e2t(ji,jj) &
& * MAX( 0.e0, rn2(ji,jj,jk) ) * zav_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
END DO
END DO
END DO
ztpc= rho0 / ( rn_tfe * rn_me ) * ztpc
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' N Total power consumption by av_tide : ztpc = ', ztpc * 1.e-12 ,'TW'
ENDIF
! ! ----------------------- !
! ! ITF tidal mixing ! (update zav_tide)
! ! ----------------------- !
IF( ln_tmx_itf ) CALL tmx_itf( kt, zav_tide , Kmm)
! ! ----------------------- !
! ! Update mixing coefs !
! ! ----------------------- !
DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing
DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
DO ji = 1, jpi
p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
END DO
END DO
END DO
! !* output tidal mixing coefficient
CALL iom_put( "av_tmx", zav_tide )
IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab3d_1=zav_tide , clinfo1=' tmx - av_tide: ', tab3d_2=p_avt, clinfo2=' p_avt: ', kdim=jpk)
!
END SUBROUTINE zdf_tmx
SUBROUTINE tmx_itf( kt, pav , Kmm)
!!----------------------------------------------------------------------
!! *** ROUTINE tmx_itf ***
!!
!! ** Purpose : modify the vertical eddy diffusivity coefficients
!! (pav) in the Indonesian Through Flow area (ITF).
!!
!! ** Method : - Following Koch-Larrouy et al. (2007), in the ITF defined
!! by msk_itf (read in a file, see tmx_init), the tidal
!! mixing coefficient is computed with :
!! * q=1 (i.e. all the tidal energy remains trapped in
!! the area and thus is used for mixing)
!! * the vertical distribution of the tifal energy is a
!! proportional to N above the thermocline (d(N^2)/dz > 0)
!! and to N^2 below the thermocline (d(N^2)/dz < 0)
!!
!! ** Action : av_tide updated in the ITF area (msk_itf)
!!
!! References : Koch-Larrouy et al. 2007, GRL
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt, Kmm ! ocean time-step
REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pav ! Tidal mixing coef.
!!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zcoef, ztpc ! temporary scalar
REAL(wp), DIMENSION(jpi,jpj) :: zkz ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zsum1 , zsum2 , zsum ! - -
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zempba_3d_1, zempba_3d_2 ! 3D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zempba_3d , zdn2dz ! - -
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zavt_itf ! - -
!!----------------------------------------------------------------------
!
! ! compute the form function using N2 at each time step
zdn2dz (:,:,jpk) = 0.e0
zempba_3d_1(:,:,jpk) = 0.e0
zempba_3d_2(:,:,jpk) = 0.e0
DO jk = 1, jpkm1
zdn2dz (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1) ! Vertical profile of dN2/dz
!CDIR NOVERRCHK
zempba_3d_1(:,:,jk) = SQRT( MAX( 0.e0, rn2(:,:,jk) ) ) ! - - of N
zempba_3d_2(:,:,jk) = MAX( 0.e0, rn2(:,:,jk) ) ! - - of N^2
END DO
!
zsum (:,:) = 0.e0
zsum1(:,:) = 0.e0
zsum2(:,:) = 0.e0
DO jk= 2, jpk
zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * e3w(:,:,jk,Kmm) * tmask(:,:,jk) * tmask(:,:,jk-1)
zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * e3w(:,:,jk,Kmm) * tmask(:,:,jk) * tmask(:,:,jk-1)
END DO
DO jj = 1, jpj
DO ji = 1, jpi
IF( zsum1(ji,jj) /= 0.e0 ) zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
IF( zsum2(ji,jj) /= 0.e0 ) zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)
END DO
END DO
DO jk= 1, jpk
DO jj = 1, jpj
DO ji = 1, jpi
zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) ) ! =0 if dN2/dz > 0, =1 otherwise
ztpc = zempba_3d_1(ji,jj,jk) * zsum1(ji,jj) * zcoef &
& + zempba_3d_2(ji,jj,jk) * zsum2(ji,jj) * ( 1. - zcoef )
!
zempba_3d(ji,jj,jk) = ztpc
zsum (ji,jj) = zsum(ji,jj) + ztpc * e3w(ji,jj,jk,Kmm)
END DO
END DO
END DO
DO jj = 1, jpj
DO ji = 1, jpi
IF( zsum(ji,jj) > 0.e0 ) zsum(ji,jj) = 1.e0 / zsum(ji,jj)
END DO
END DO
! ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min)
zcoef = rn_tfe_itf / ( rn_tfe * rho0 )
DO jk = 1, jpk
zavt_itf(:,:,jk) = MIN( 10.e-4, zcoef * en_tmx(:,:) * zsum(:,:) * zempba_3d(:,:,jk) &
& / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk) )
END DO
zkz(:,:) = 0.e0 ! Associated potential energy consummed over the whole water column
DO jk = 2, jpkm1
zkz(:,:) = zkz(:,:) + e3w(:,:,jk,Kmm) * MAX( 0.e0, rn2(:,:,jk) ) * rho0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
END DO
DO jj = 1, jpj ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
DO ji = 1, jpi
IF( zkz(ji,jj) /= 0.e0 ) zkz(ji,jj) = en_tmx(ji,jj) * rn_tfe_itf / rn_tfe / zkz(ji,jj)
END DO
END DO
DO jk = 2, jpkm1 ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1) ! kz max = 120 cm2/s
END DO
IF( kt == nit000 ) THEN ! diagnose the nergy consumed by zavt_itf
ztpc = 0.e0
DO jk= 1, jpk
DO jj= 1, jpj
DO ji= 1, jpi
ztpc = ztpc + e1t(ji,jj) * e2t(ji,jj) * e3w(ji,jj,jk,Kmm) * MAX( 0.e0, rn2(ji,jj,jk) ) &
& * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
END DO
END DO
END DO
ztpc= rho0 * ztpc / ( rn_me * rn_tfe_itf )
IF(lwp) WRITE(numout,*) ' N Total power consumption by zavt_itf: ztpc = ', ztpc * 1.e-12 ,'TW'
ENDIF
! ! Update pav with the ITF mixing coefficient
DO jk = 2, jpkm1
pav(:,:,jk) = pav (:,:,jk) * ( 1.e0 - mask_itf(:,:) ) &
& + zavt_itf(:,:,jk) * mask_itf(:,:)
END DO
!
END SUBROUTINE tmx_itf
SUBROUTINE zdf_tmx_init
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_tmx_init ***
!!
!! ** Purpose : Initialization of the vertical tidal mixing, Reading
!! of M2 and K1 tidal energy in nc files
!!
!! ** Method : - Read the namtmx namelist and check the parameters
!!
!! - Read the input data in NetCDF files :
!! M2 and K1 tidal energy. The total tidal energy, en_tmx,
!! is the sum of M2, K1 and S2 energy where S2 is assumed
!! to be: S2=(1/2)^2 * M2
!! mask_itf, a mask array that determine where substituing
!! the standard Simmons et al. (2005) formulation with the
!! one of Koch_Larrouy et al. (2007).
!!
!! - Compute az_tmx, a 3D coefficient that allows to compute
!! the standard tidal-induced vertical mixing as follows:
!! Kz_tides = az_tmx / max( rn_n2min, N^2 )
!! with az_tmx a bottom intensified coefficient is given by:
!! az_tmx(z) = en_tmx / ( rho0 * rn_htmx ) * EXP( -(H-z)/rn_htmx )
!! / ( 1. - EXP( - H /rn_htmx ) )
!! where rn_htmx the characteristic length scale of the bottom
!! intensification, en_tmx the tidal energy, and H the ocean depth
!!
!! ** input : - Namlist namtmx
!! - NetCDF file : M2_ORCA2.nc, K1_ORCA2.nc, and mask_itf.nc
!!
!! ** Action : - Increase by 1 the nstop flag is setting problem encounter
!! - defined az_tmx used to compute tidal-induced mixing
!!
!! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
!! Koch-Larrouy et al. 2007, GRL.
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: inum ! local integer
INTEGER :: ios
REAL(wp) :: ztpc, ze_z ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: zem2, zek1 ! read M2 and K1 tidal energy
REAL(wp), DIMENSION(jpi,jpj) :: zkz ! total M2, K1 and S2 tidal energy
REAL(wp), DIMENSION(jpi,jpj) :: zfact ! used for vertical structure function
REAL(wp), DIMENSION(jpi,jpj) :: zhdep ! Ocean depth
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpc ! power consumption
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zav_tide ! tidal mixing coefficient
!!
NAMELIST/namzdf_tmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf
!!----------------------------------------------------------------------
!
! Namelist namzdf_tmx in reference namelist : Tidal Mixing
READ ( numnam_ref, namzdf_tmx, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist' )
! Namelist namzdf_tmx in configuration namelist : Tidal Mixing
READ ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist' )
IF(lwm) WRITE ( numond, namzdf_tmx )
IF(lwp) THEN ! Control print
WRITE(numout,*)
WRITE(numout,*) 'zdf_tmx_init : tidal mixing'
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namzdf_tmx : set tidal mixing parameters'
WRITE(numout,*) ' Vertical decay scale for turbulence = ', rn_htmx
WRITE(numout,*) ' Brunt-Vaisala frequency threshold = ', rn_n2min
WRITE(numout,*) ' Tidal dissipation efficiency = ', rn_tfe
WRITE(numout,*) ' Mixing efficiency = ', rn_me
WRITE(numout,*) ' ITF specific parameterisation = ', ln_tmx_itf
WRITE(numout,*) ' ITF tidal dissipation efficiency = ', rn_tfe_itf
ENDIF
! ! allocate tmx arrays
IF( zdf_tmx_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
IF( ln_tmx_itf ) THEN ! read the Indonesian Through Flow mask
CALL iom_open('mask_itf',inum)
CALL iom_get (inum, jpdom_global, 'tmaskitf',mask_itf,1) !
CALL iom_close(inum)
ENDIF
! read M2 tidal energy flux : W/m2 ( zem2 < 0 )
CALL iom_open('M2rowdrg',inum)
CALL iom_get (inum, jpdom_global, 'field',zem2,1) !
CALL iom_close(inum)
! read K1 tidal energy flux : W/m2 ( zek1 < 0 )
CALL iom_open('K1rowdrg',inum)
CALL iom_get (inum, jpdom_global, 'field',zek1,1) !
CALL iom_close(inum)
! Total tidal energy ( M2, S2 and K1 with S2=(1/2)^2 * M2 )
! only the energy available for mixing is taken into account,
! (mixing efficiency tidal dissipation efficiency)
en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * ssmask(:,:)
!============
!TG: Bug for VVL? Should this section be moved out of _init and be updated at every timestep?
! Vertical structure (az_tmx)
DO jj = 1, jpj ! part independent of the level
DO ji = 1, jpi
zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean
zfact(ji,jj) = rho0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
END DO
END DO
DO jk= 1, jpk ! complete with the level-dependent part
DO jj = 1, jpj
DO ji = 1, jpi
az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-gdepw_0(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
END DO
END DO
END DO
!===========
IF( lwp ) THEN
! Control print
! Total power consumption due to vertical mixing
! zpc = rho0 * 1/rn_me * rn2 * zav_tide
zav_tide(:,:,:) = 0.e0
DO jk = 2, jpkm1
zav_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
END DO
ztpc = 0.e0
zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:)
DO jk= 2, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
ztpc = ztpc + e3w_0(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
END DO
END DO
END DO
ztpc= rho0 * 1/(rn_tfe * rn_me) * ztpc
WRITE(numout,*)
WRITE(numout,*) ' Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
! control print 2
zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )
zkz(:,:) = 0.e0
DO jk = 2, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
zkz(ji,jj) = zkz(ji,jj) + e3w_0(ji,jj,jk) * MAX(0.e0, rn2(ji,jj,jk)) * rho0 * zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
END DO
END DO
END DO
! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
DO jj = 1, jpj
DO ji = 1, jpi
IF( zkz(ji,jj) /= 0.e0 ) THEN
zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
ENDIF
END DO
END DO
ztpc = 1.e50
DO jj = 1, jpj
DO ji = 1, jpi
IF( zkz(ji,jj) /= 0.e0 ) THEN
ztpc = Min( zkz(ji,jj), ztpc)
ENDIF
END DO
END DO
WRITE(numout,*) ' Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) )
DO jk = 2, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s
END DO
END DO
END DO
ztpc = 0.e0
zpc(:,:,:) = Max(0.e0,rn2(:,:,:)) * zav_tide(:,:,:)
DO jk= 1, jpk
DO jj = 1, jpj
DO ji = 1, jpi
ztpc = ztpc + e3w_0(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
END DO
END DO
END DO
ztpc= rho0 * 1/(rn_tfe * rn_me) * ztpc
WRITE(numout,*) ' 2 Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
DO jk = 1, jpk
ze_z = SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk) * tmask_i(:,:) ) &
& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
ztpc = 1.E50
DO jj = 1, jpj
DO ji = 1, jpi
IF( zav_tide(ji,jj,jk) /= 0.e0 ) ztpc =Min( ztpc, zav_tide(ji,jj,jk) )
END DO
END DO
WRITE(numout,*) ' N2 min - jk= ', jk,' ', ze_z * 1.e4,' cm2/s min= ',ztpc*1.e4, &
& 'max= ', MAXVAL(zav_tide(:,:,jk) )*1.e4, ' cm2/s'
END DO
WRITE(numout,*) ' e_tide : ', SUM( e1t*e2t*en_tmx ) / ( rn_tfe * rn_me ) * 1.e-12, 'TW'
WRITE(numout,*)
WRITE(numout,*) ' Initial profile of tidal vertical mixing'
DO jk = 1, jpk
DO jj = 1,jpj
DO ji = 1,jpi
zkz(ji,jj) = az_tmx(ji,jj,jk) /MAX( rn_n2min, rn2(ji,jj,jk) )
END DO
END DO
ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) &
& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
WRITE(numout,*) ' jk= ', jk,' ', ze_z * 1.e4,' cm2/s'
END DO
DO jk = 1, jpk
zkz(:,:) = az_tmx(:,:,jk) /rn_n2min
ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) &
& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
WRITE(numout,*)
WRITE(numout,*) ' N2 min - jk= ', jk,' ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4, &
& 'max= ', MAXVAL(zkz)*1.e4, ' cm2/s'
END DO
!
ENDIF
!
END SUBROUTINE zdf_tmx_init
!!======================================================================
END MODULE zdftmx
src/OCE/nemogcm.F90
View file @ 4a9f4b18
......@@ -73,6 +73,7 @@ MODULE nemogcm
USE lib_mpp ! distributed memory computing
USE mppini ! shared/distributed memory setting (mpp_init routine)
USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
USE sbccpl
USE halo_mng ! halo manager
IMPLICIT NONE
......@@ -173,6 +174,9 @@ CONTAINS
IF ( istp == nitend ) elapsed_time = zstptiming - elapsed_time
ENDIF
!
IF (lk_oasis) THEN
CALL sbc_cpl_snd( istp, Nbb, Nnn ) ! Coupling to atmos
ENDIF
# if defined key_qco || defined key_linssh
CALL stp_MLF( istp )
# else
......@@ -262,7 +266,7 @@ CONTAINS
#if defined key_xios
IF( Agrif_Root() ) THEN
IF( lk_oasis ) THEN
CALL cpl_init( "oceanx", ilocal_comm ) ! nemo local communicator given by oasis
CALL cpl_init( "toyoce", ilocal_comm ) ! nemo local communicator given by oasis
CALL xios_initialize( "not used" , local_comm =ilocal_comm ) ! send nemo communicator to xios
ELSE
CALL xios_initialize( "for_xios_mpi_id", return_comm=ilocal_comm ) ! nemo local communicator given by xios
......@@ -272,7 +276,7 @@ CONTAINS
#else
IF( lk_oasis ) THEN
IF( Agrif_Root() ) THEN
CALL cpl_init( "oceanx", ilocal_comm ) ! nemo local communicator given by oasis
CALL cpl_init( "toyoce", ilocal_comm ) ! nemo local communicator given by oasis
ENDIF
CALL mpp_start( ilocal_comm )
ELSE
......@@ -496,6 +500,11 @@ CONTAINS
!
IF(lwp) WRITE(numout,cform_aaa) ! Flag AAAAAAA
!
IF (nstop > 0) THEN
CALL CTL_STOP('STOP','Critical errors in NEMO initialisation')
END IF
IF( ln_timing ) CALL timing_stop( 'nemo_init')
!
END SUBROUTINE nemo_init
......
src/OCE/par_oce.F90
View file @ 4a9f4b18
......@@ -99,6 +99,12 @@ MODULE par_oce
INTEGER, PUBLIC :: Ni_0, Nj_0 !: local domain size without halo
INTEGER, PUBLIC :: Ni0glo, Nj0glo !: global domain size without halo
INTEGER, PUBLIC :: Nis0_ext !: start I-index with wrap/N-fold
INTEGER, PUBLIC :: Nie0_ext !: end I-index with wrap/N-fold
INTEGER, PUBLIC :: Njs0_ext !: start J-index with wrap/N-fold
INTEGER, PUBLIC :: Nje0_ext !: end J-index with wrap/N-fold
INTEGER, PUBLIC :: Ni_0_ext, Nj_0_ext !: local domain size with wrap/N-fold
INTEGER, PUBLIC :: Ni0glo_ext, Nj0glo_ext !: global domain size with wrap/N-fold
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: par_oce.F90 15119 2021-07-13 14:43:22Z jchanut $
......
src/OCE/step.F90
View file @ 4a9f4b18
......@@ -172,6 +172,7 @@ CONTAINS
CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency
! VERTICAL PHYSICS
! lbc_lnk needed for zdf_sh2 when using nn_hls = 2, moved here to allow tiling in zdf_phy
IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp', avm_k, 'W', 1.0_wp )
......@@ -296,6 +297,7 @@ CONTAINS
IF( ln_diadct ) CALL dia_dct ( kstp, Nnn ) ! Transports
CALL dia_ar5 ( kstp, Nnn ) ! ar5 diag
CALL dia_ptr ( kstp, Nnn ) ! Poleward adv/ldf TRansports diagnostics
CALL dia_prod ( kstp, Nnn ) ! ocean model: products
CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs
IF( ln_crs ) CALL crs_fld ( kstp, Nnn ) ! ocean model: online field coarsening & output
IF( lk_diadetide ) CALL dia_detide( kstp ) ! Weights computation for daily detiding of model diagnostics
......@@ -429,7 +431,7 @@ CONTAINS
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Coupled mode
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
!!IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
!
#if defined key_xios
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
......
src/OCE/step_oce.F90
View file @ 4a9f4b18
......@@ -85,6 +85,7 @@ MODULE step_oce
USE diahsb ! heat, salt and volume budgets (dia_hsb routine)
USE diacfl ! CFL diagnostics (dia_cfl routine)
USE diaobs ! Observation operator (dia_obs routine)
USE diaprod
USE diadetide ! Weights computation for daily detiding of model diagnostics
USE diamlr ! IOM context management for multiple-linear-regression analysis
USE flo_oce ! floats variables
......