diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index 5c329d3aca36475e3b32fe3224e2aa7ccdc91705..e52985c3c130d84046925d3ce20ec87dfab98cde 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -323,6 +323,7 @@
ln_rstart_abl = .false.
ln_hpgls_frc = .false.
+ ln_pga_abl = .false. ! ABL pressure gradient anomaly forcing
ln_geos_winds = .false.
ln_smth_pblh = .false.
nn_dyn_restore = 0 ! restoring option for dynamical ABL variables: = 0 no restoring
diff --git a/src/ABL/ablmod.F90 b/src/ABL/ablmod.F90
index d14691e4e56216314d10e8143d85aca8b7cace4c..633e24652adcf1925972bed0d3c0cb0049455860 100644
--- a/src/ABL/ablmod.F90
+++ b/src/ABL/ablmod.F90
@@ -105,6 +105,8 @@ CONTAINS
REAL(wp), DIMENSION(A2D(0)) :: zwnd_i, zwnd_j
REAL(wp), DIMENSION(A2D(0)) :: zsspt, ztabs, zpre
REAL(wp), DIMENSION(A1Di(0),2:jpka) :: zCF
+ REAL(wp), DIMENSION(A2D(1),1:jpka) :: zpdif
+ REAL(wp), DIMENSION(A2D(0),1:jpka) :: zpabl, zpdta, zpgau, zpgav
REAL(wp), DIMENSION(A1Di(0),1:jpka) :: z_elem_a, z_elem_b, z_elem_c
!
INTEGER :: ji, jj, jk, jtra, jbak ! dummy loop indices
@@ -140,10 +142,6 @@ CONTAINS
!zrough(:,:) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(:,:), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce
zrough(:,:) = z0_from_Cd( ght_abl(2), pCd_du(:,:) / MAX( pwndm(:,:), 0.5_wp ) ) ! #LB: z0_from_Cd is define in sbc_phy.F90...
- ! sea surface potential temperature [K]
- !zsspt(:,:) = theta_exner( psst(A2D(0))+rt0, pslp_dta(:,:) )
-
-
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! 1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -282,11 +280,39 @@ CONTAINS
END IF
!
IF( ln_hpgls_frc ) THEN
+
+ IF( ln_pga_abl ) THEN
+
+ DO_3D( 0, 0, 0, 0, 0, jpka)
+ zpabl(ji,jj,jk) = pres_temp( tq_abl(ji,jj,jk,nt_n,jp_qa), pslp_dta(ji,jj), ght_abl(jk), ptpot=tq_abl(ji,jj,jk,nt_n,jp_ta) )
+ zpdta(ji,jj,jk) = pres_temp( pq_dta(ji,jj,jk) , pslp_dta(ji,jj), ght_abl(jk), ptpot=pt_dta(ji,jj,jk) )
+ zpdif(ji,jj,jk) = zpabl(ji,jj,jk) - zpdta(ji,jj,jk)
+ END_3D
+ zpdif(:,:,:) = MAX( MIN( zpdif(:,:,:), 10._wp ), -10._wp ) ! limiter due to ice/oce inconsistencies between nemo/si3 & atm frc
+
+ DO jk = 1, jpka
+ CALL smooth_pblh( zpdif(:,:,:), msk_abl(:,:) )
+ END DO
+ CALL lbc_lnk( 'ablmod', zpdif(:,:,:), 'T', 1.0_wp)
+
+ DO_3D( 0, 0, 0, 0, 1, jpka)
+ zpgau(ji,jj,jk) = ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) * &
+ & ( (zpdif(ji+1,jj,jk) - zpdif(ji ,jj,jk)) * r1_e1u(ji ,jj) * umask(ji ,jj,1) + &
+ & (zpdif(ji ,jj,jk) - zpdif(ji-1,jj,jk)) * r1_e1u(ji-1,jj) * umask(ji-1,jj,1) )
+ zpgav(ji,jj,jk) = ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) * &
+ & ( (zpdif(ji,jj+1,jk) - zpdif(ji,jj ,jk)) * r1_e2v(ji,jj ) * vmask(ji,jj ,1) + &
+ & (zpdif(ji,jj ,jk) - zpdif(ji,jj-1,jk)) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,1) )
+ END_3D
+ ELSE
+ zpgau(:,:,:) = 0._wp
+ zpgav(:,:,:) = 0._wp
+ ENDIF
+
DO_1Dj( 0, 0 ) ! outer loop
DO jk = 1, jpka
DO_1Di( 0, 0 )
- u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk)
- v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)
+ u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * ( pgu_dta(ji,jj,jk) + zpgau(ji,jj,jk) )
+ v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * ( pgv_dta(ji,jj,jk) + zpgav(ji,jj,jk) )
END_1D
ENDDO
END_1D
@@ -788,8 +814,9 @@ CONTAINS
END DO
!!FL should not be needed because of Patankar procedure
- tke_abl(A2D(0),1:jpka,nt_a) = MAX( tke_abl(A2D(0),1:jpka,nt_a), tke_min )
-
+ DO_1Di( 0, 0 )
+ tke_abl(ji,jj,1:jpka,nt_a) = MAX( tke_abl(ji,jj,1:jpka,nt_a), tke_min )
+ END_1D
!!
!! Diagnose PBL height
!! ----------------------------------------------------------
diff --git a/src/ABL/par_abl.F90 b/src/ABL/par_abl.F90
index 35dcda6f381faf91b4e6bd1ba9a7db8020dde3e5..88f1995b833c4a827e94acabdd738311a13ba686 100644
--- a/src/ABL/par_abl.F90
+++ b/src/ABL/par_abl.F90
@@ -26,6 +26,7 @@ MODULE par_abl
INTEGER , PUBLIC :: nn_dyn_restore !: restoring option for dynamical ABL variables
LOGICAL , PUBLIC :: ln_geos_winds !: large-scale restoring of ABL winds toward geostrophic winds
LOGICAL , PUBLIC :: ln_hpgls_frc !: forcing of ABL winds by large-scale pressure gradient
+ LOGICAL , PUBLIC :: ln_pga_abl !: ABL pressure gradient anomaly forcing
LOGICAL , PUBLIC :: ln_smth_pblh !: smoothing of atmospheric PBL height
!LOGICAL , PUBLIC :: ln_topbc_neumann = .FALSE. !: idealised testcases only
diff --git a/src/ABL/sbcabl.F90 b/src/ABL/sbcabl.F90
index f4ac4657b3e061b336fdaaa2f35a30c451ef4bc9..d2efc748e104bdcd9145a54475a033475aa46153 100644
--- a/src/ABL/sbcabl.F90
+++ b/src/ABL/sbcabl.F90
@@ -73,7 +73,7 @@ CONTAINS
& ln_hpgls_frc, ln_geos_winds, nn_dyn_restore, &
& rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max, &
& nn_amxl, rn_Cm, rn_Ct, rn_Ce, rn_Ceps, rn_Rod, rn_Ric, &
- & rn_vfac, ln_smth_pblh
+ & rn_vfac, ln_smth_pblh, ln_pga_abl
!!---------------------------------------------------------------------
! Namelist namsbc_abl in reference namelist : ABL parameters
@@ -107,6 +107,9 @@ CONTAINS
ln_geos_winds = .FALSE.
WRITE(numout,*) ' ABL -- geostrophic guide disabled (not compatible with ln_hpgls_frc = .T.)'
END IF
+ IF( ln_pga_abl ) THEN
+ WRITE(numout,*) ' ABL -- pressure gradient anomaly forcing'
+ END IF
ELSE IF( ln_geos_winds ) THEN
WRITE(numout,*) ' ABL -- winds forced by geostrophic winds'
ELSE
@@ -272,11 +275,6 @@ CONTAINS
nt_n = 1; nt_a = 2
! initialize ABL from data or restart
- u_abl (:,:,:,nt_a ) = 0._wp
- v_abl (:,:,:,nt_a ) = 0._wp
- tq_abl (:,:,:,nt_a,: ) = 0._wp
- tke_abl(:,:,:,nt_a ) = 0._wp
-
IF( ln_rstart_abl ) THEN
CALL abl_rst_read
ELSE
@@ -320,11 +318,9 @@ CONTAINS
!!---------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! ocean time step
!!
- !REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp
- REAL(wp), DIMENSION(A2D(0)) :: zssq, zcd_du, zsen, zlat, zevp
+ REAL(wp), DIMENSION(A2D(0)) :: zssq, zcd_du, zsen, zlat, zevp
#if defined key_si3
- !REAL(wp), DIMENSION(jpi,jpj) :: zssqi, zcd_dui, zseni, zevpi
- REAL(wp), DIMENSION(A2D(0)) :: zssqi, zcd_dui, zseni, zevpi
+ REAL(wp), DIMENSION(A2D(0)) :: zssqi, zcd_dui, zseni, zevpi
#endif
INTEGER :: jbak, jbak_dta, ji, jj
!!---------------------------------------------------------------------