Skip to content
Snippets Groups Projects
Commit c47e4fd4 authored by Richard Hill's avatar Richard Hill
Browse files

Various fixes to NEMO applicable to forced ocean/seaice and coupled models

parent 9953c80b
No related branches found
No related tags found
No related merge requests found
......@@ -453,6 +453,11 @@ MODULE ice
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: qcn_ice_top !: Surface conduction flux (W/m2)
!
!!----------------------------------------------------------------------
!! * Only for atmospheric coupling
!!----------------------------------------------------------------------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Ice fractional area at last coupling time
!
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: ice.F90 15388 2021-10-17 11:33:47Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
......@@ -551,6 +556,10 @@ CONTAINS
ii = ii + 1
ALLOCATE( t_si(jpi,jpj,jpl) , tm_si(jpi,jpj) , qcn_ice_bot(jpi,jpj,jpl) , qcn_ice_top(jpi,jpj,jpl) , STAT = ierr(ii) )
! * For atmospheric coupling
ii = ii + 1
ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(ii) )
ice_alloc = MAXVAL( ierr(:) )
IF( ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' )
!
......
......@@ -996,24 +996,22 @@ CONTAINS
!!-----------------------------------------------------------------------
!! *** ROUTINE ice_dyn_1d2d ***
!!
!! ** Purpose : move arrays from 1d to 2d and the reverse
!! ** Purpose : move arrays between 1d <=> 2d forms and
!! between 2d <=> 3d forms
!!-----------------------------------------------------------------------
INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
INTEGER, INTENT(in) :: kn ! 1= from 2D/3D to 1D/2D
! 2= from 1D/2D to 2D/3D
!
INTEGER :: jl, jk ! dummy loop indices
!!-----------------------------------------------------------------------
!
SELECT CASE( kn )
! !---------------------!
CASE( 1 ) !== from 2D to 1D ==!
! !---------------------!
! !---------------------------!
CASE( 1 ) !== from 2D/3D to 1D/2D ==!
! !---------------------------!
! fields used but not modified
CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
! the following fields are modified in this routine
!!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
!!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
!!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
......@@ -1035,9 +1033,9 @@ CONTAINS
CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) )
!
! !---------------------!
CASE( 2 ) !== from 1D to 2D ==!
! !---------------------!
! !---------------------------!
CASE( 2 ) !== from 1D/2D to 2D/3D ==!
! !---------------------------!
CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
......
......@@ -324,6 +324,11 @@ CONTAINS
ENDIF
ENDIF
! If this is a coupled model we need to pick up a_i for use as a_i_last_couple
IF (ln_cpl) then
a_i_last_couple = a_i
ENDIF
IF(.NOT.lrxios) CALL iom_delay_rst( 'READ', 'ICE', numrir ) ! read only ice delayed global communication variables
! ! ---------------------------------- !
ELSE ! == case of a simplified restart == !
......
......@@ -32,14 +32,15 @@ MODULE icetab
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab1d, tab2d )
SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab2d, tab3d )
!!----------------------------------------------------------------------
!! *** ROUTINE tab_2d_1d ***
!! *** ROUTINE tab_3d_2d ***
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ndim1d ! 1d size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in ) :: tab2d ! input 2D field
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) :: tab1d ! output 1D field
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in ) :: tab3d ! input 3D field
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) :: tab2d ! output 2D field
!
INTEGER :: jl, jn, jid, jjd
!!----------------------------------------------------------------------
......@@ -47,7 +48,7 @@ CONTAINS
DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1
jjd = ( tab_ind(jn) - 1 ) / jpi + 1
tab1d(jn,jl) = tab2d(jid,jjd,jl)
tab2d(jn,jl) = tab3d(jid,jjd,jl)
END DO
END DO
END SUBROUTINE tab_3d_2d
......@@ -72,14 +73,14 @@ CONTAINS
END SUBROUTINE tab_2d_1d
SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab1d, tab2d )
SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab2d, tab3d )
!!----------------------------------------------------------------------
!! *** ROUTINE tab_2d_1d ***
!! *** ROUTINE tab_2d_3d ***
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ndim1d ! 1D size
INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab1d ! input 1D field
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) :: tab2d ! output 2D field
REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab2d ! input 2D field
REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) :: tab3d ! output 3D field
!
INTEGER :: jl, jn, jid, jjd
!!----------------------------------------------------------------------
......@@ -87,7 +88,7 @@ CONTAINS
DO jn = 1, ndim1d
jid = MOD( tab_ind(jn) - 1 , jpi ) + 1
jjd = ( tab_ind(jn) - 1 ) / jpi + 1
tab2d(jid,jjd,jl) = tab1d(jn,jl)
tab3d(jid,jjd,jl) = tab2d(jn,jl)
END DO
END DO
END SUBROUTINE tab_2d_3d
......
......@@ -692,7 +692,7 @@ CONTAINS
!! - vertical interpolation: simple averaging
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) :: pe3_out ! output interpolated e3
CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors
! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
!
......
......@@ -167,7 +167,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
REAL(wp), DIMENSION(jpi,jpj) :: zhV! fluxes
REAL(dp), DIMENSION(jpi,jpj) :: zhU! fluxes
REAL(wp), DIMENSION(jpi,jpj) :: zhU! fluxes
!!st#if defined key_qco
!!st REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
!!st#endif
......@@ -270,6 +270,11 @@ CONTAINS
zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
END_2D
ELSE
! Ensure zhU and zhV are initialised to SOMETHING at all points to avoid referencing
! uninitialsed values in halos later on!
zhU(:,:) = 0._wp
zhV(:,:) = 0._wp
ENDIF
!
! != Add bottom stress contribution from baroclinic velocities =!
......@@ -502,13 +507,13 @@ CONTAINS
IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
!
! ! resulting flux at mid-step (not over the full domain)
DO_2D( 1, 0, 1, 1 ) ! not jpi-column
zhU(ji,jj) = e2u(ji,jj) * ua_e(ji,jj) * zhup2_e(ji,jj)
END_2D
DO_2D( 1, 1, 1, 0 ) ! not jpj-row
zhV(ji,jj) = e1v(ji,jj) * va_e(ji,jj) * zhvp2_e(ji,jj)
END_2D
!
#if defined key_agrif
! Set fluxes during predictor step to ensure volume conservation
IF( ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
......@@ -519,6 +524,12 @@ CONTAINS
CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
!
ENDIF
! It seems safest to do this here since zhU and zhV are not initially calculated in halos
! by this code or by wad_Umsk, but halo values (ji-1 and jj-1) ARE required in the zhdiv
! sea level calculation. The current trunk (Feb 2024) has resolved all these issues by rewriting.
CALL lbc_lnk( 'dynspg_ts', zhU, 'U', -1._wp)
CALL lbc_lnk( 'dynspg_ts', zhV, 'V', -1._wp)
!
!
! Compute Sea Level at step jit+1
......@@ -529,15 +540,16 @@ CONTAINS
zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( ssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
END_2D
!
CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._dp, zhU, 'U', -1._dp)
CALL lbc_lnk( 'dynspg_ts', zhV, 'V', -1._wp )
CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._dp)
!
! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
IF( ln_bdy ) CALL bdy_ssh( ssha_e )
#if defined key_agrif
CALL agrif_ssh_ts( jn )
#endif
!
! ! Sum over sub-time-steps to compute advective velocities
za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
......@@ -556,6 +568,7 @@ CONTAINS
!
! Sea Surface Height at u-,v-points (vvl case only)
IF( .NOT.ln_linssh ) THEN
#if defined key_qcoTest_FluxForm
! ! 'key_qcoTest_FluxForm' : simple ssh average
DO_2D( 1, 0, 1, 1 )
......@@ -623,6 +636,7 @@ CONTAINS
!-- u = m+1 | h * u + delta_t' * \ h * (1-r)*g * grad_x( ssh') - h * f * k vect u + h * frc / | --!
!-- h \ / --!
!------------------------------------------------------------------------------------------------------------------------!
IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
DO_2D( 0, 0, 0, 0 )
ua_e(ji,jj) = ( un_e(ji,jj) &
......@@ -676,7 +690,7 @@ CONTAINS
ENDIF
IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
DO_2D( 0, 0, 0, 0 )
DO_2D( 0, 0, 0, 0 )
hu_e (ji,jj) = hu_0(ji,jj) + zsshu_a(ji,jj)
hur_e(ji,jj) = ssumask(ji,jj) / ( hu_e(ji,jj) + 1._wp - ssumask(ji,jj) )
hv_e (ji,jj) = hv_0(ji,jj) + zsshv_a(ji,jj)
......@@ -686,10 +700,10 @@ CONTAINS
!
IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
, hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
, hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
ELSE
CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
ENDIF
! ! open boundaries
IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
......@@ -728,7 +742,6 @@ CONTAINS
! ! Sum sea level
pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
! ! ==================== !
END DO ! end loop !
! ! ==================== !
! -----------------------------------------------------------------------------
......@@ -1085,7 +1098,13 @@ CONTAINS
! ! Allocate time-splitting arrays
IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
!
! ! read restart when needed
! RSRH. I've just copied the following from the current NEMO main. It shouldnt affect evolution (!) but
! does help with debugging and testing!
! init some arrays for debug sette
ssha_e(:,:) = 0._wp
!
! !: restart/initialise
CALL ts_rst( nit000, 'READ' )
!
END SUBROUTINE dyn_spg_ts_init
......@@ -1165,8 +1184,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ji ,jj ! dummy loop indices
REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhV
REAL(dp), DIMENSION(jpi,jpj), INTENT(in ) :: zhU
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
!!----------------------------------------------------------------------
SELECT CASE( nvor_scheme )
......@@ -1272,8 +1290,7 @@ CONTAINS
!! ** Action : ptmsk : wetting & drying t-mask
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phV, pu, pv! ocean velocities and transports
REAL(dp), DIMENSION(jpi,jpj), INTENT(inout) :: phU! ocean velocities and transports
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
!
INTEGER :: ji, jj ! dummy loop indices
......
......@@ -1436,7 +1436,7 @@ ENDIF
! This is rather confused by the fact that jpjglo has a value BIGGER
! than it did at pre 4.2... e.g. for ORCA1 it's set to 333 instead of 332
! which is rather baffling and confuses some dimensioning calculations
! if we're not very very careful
! which previously worked fine pre 4.2.
! Set up dimensions for old style coupling exchanges on extended grid
......
......@@ -45,7 +45,6 @@ MODULE cpl_oasis3
PRIVATE
#if ! defined key_oasis3
! Dummy interface to oasis_get if not using oasis
! RSRH Is this really needed
INTERFACE oasis_get
MODULE PROCEDURE oasis_get_1d, oasis_get_2d
END INTERFACE
......@@ -253,8 +252,8 @@ CONTAINS
CALL oasis_def_partition ( id_part_2d, paral, nerror, Ni0glo*Nj0glo ) ! global number of points, excluding halos
! RSRH Set up 2D box partition for compatibility with existing weights files
! so we DONT HAVE TO GENERATE AND MANAGE MULTIPLE SETS OF WEIGHTS PURLEY BECAUSE OF AN
! ARBITRARY CHANGE IN THE SOURCE CODE!
! so we don't have to generate and manage multiple sets of weights purely because of
! the changes to nemo 4.2+ code!
! This is just a hack for global cyclic models for the time being
Ni0glo_ext = jpiglo
......@@ -569,13 +568,14 @@ CONTAINS
INTEGER , INTENT( out) :: kinfo ! OASIS3 info argument
!!
INTEGER :: jc,jm ! local loop index
LOGICAL :: llaction, ll_1st
LOGICAL :: llaction, ll_1st, lrcv
!!--------------------------------------------------------------------
!
! receive local data from OASIS3 on every process
!
kinfo = OASIS_idle
lrcv=.FALSE.
!
DO jc = 1, srcv(kid)%nct
......@@ -594,9 +594,10 @@ CONTAINS
& WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
IF( llaction ) THEN ! data received from oasis do not include halos
! RSRH, but DO cater for wrap columns when using pre vn 4.2 format remapping weights.
! but DO still cater for wrap columns when using pre vn4.2 compatible remapping weights.
kinfo = OASIS_Rcv
lrcv=.TRUE.
IF( ll_1st ) THEN
pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) = exfld_ext(:,:) * pmask(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jm)
......@@ -626,14 +627,17 @@ CONTAINS
ENDDO
!--- Call lbc_lnk to populate halos of received fields.
IF( .NOT. ll_1st ) THEN
CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
ENDIF
ENDDO
! RSRH I've changed this since:
! 1) it seems multi cat fields may not be properly updated in halos when called on a per
! category basis(?)
! 2) it's more efficient to have a single call (and simpler to understand) to update ALL
! categories at the same time!
!--- Call lbc_lnk to populate halos of received fields.
IF (lrcv) then
CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,:), srcv(kid)%clgrid, srcv(kid)%nsgn )
endif
!
END SUBROUTINE cpl_rcv
......
......@@ -33,7 +33,7 @@ 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
......@@ -254,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
......@@ -1460,6 +1458,14 @@ CONTAINS
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) &
......@@ -1574,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) !
! ! ================== !
......@@ -1871,6 +1883,7 @@ CONTAINS
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
......@@ -2070,13 +2083,7 @@ 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(:,:)
......@@ -2276,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)
......@@ -2702,8 +2693,12 @@ CONTAINS
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
......@@ -3069,13 +3064,20 @@ CONTAINS
*(zoty1(ji,jj)+zoty1(ji+1,jj) &
+zoty1(ji,jj-1)+zoty1(ji+1,jj-1))
END_2D
! Ensure any N fold and wrap columns are updated
CALL lbc_lnk('zotx1', ztmp1, 'V', -1.0)
CALL lbc_lnk('zoty1', ztmp2, 'U', -1.0)
ikchoix = -1
CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,zotx1,zoty1,ikchoix)
! 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
!
......
......@@ -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) )
......
......@@ -174,12 +174,12 @@ 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
IF (lk_oasis) THEN
CALL sbc_cpl_snd( istp, Nbb, Nnn ) ! Coupling to atmos
ENDIF
CALL stp ( istp )
# endif
istp = istp + 1
......
......@@ -443,7 +443,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
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
......@@ -561,6 +561,7 @@ CONTAINS
CALL lbc_lnk( 'finalize_lbc', r3u(:,:,Kaa), 'U', 1._wp, r3v(:,:,Kaa), 'V', 1._wp, &
& r3u_f(:,:), 'U', 1._wp, r3v_f(:,:), 'V', 1._wp )
ENDIF
! !* BDY open boundaries
IF( ln_bdy ) THEN
CALL bdy_tra( kt, Kbb, pts, Kaa )
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment