From 3554e17f9cfd7cdb913b2b253b50670296208e8c Mon Sep 17 00:00:00 2001
From: Guillaume Samson <guillaume.samson@mercator-ocean.fr>
Date: Wed, 8 Dec 2021 14:29:53 +0000
Subject: [PATCH] new expression to compute air and surface potential
temperatures in bulk formulae (old ticket #2632)
---
cfgs/SHARED/namelist_ref | 1 +
sette/.gitignore | 1 +
sette/{param.cfg => param.default} | 0
src/ABL/ablmod.F90 | 106 ++++---
src/ABL/par_abl.F90 | 1 +
src/ABL/sbcabl.F90 | 11 +-
src/OCE/SBC/sbc_phy.F90 | 441 ++++++++++++++++++---------
src/OCE/SBC/sbcblk.F90 | 235 +++++++-------
src/OCE/SBC/sbcblk_algo_coare3p0.F90 | 12 +-
src/OCE/SBC/sbcblk_algo_coare3p6.F90 | 12 +-
src/OCE/SBC/sbcblk_algo_ecmwf.F90 | 11 +-
src/OCE/SBC/sbccpl.F90 | 8 +-
12 files changed, 512 insertions(+), 327 deletions(-)
rename sette/{param.cfg => param.default} (100%)
diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index a397a232e..4fa3957cb 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -330,6 +330,7 @@
rn_ldyn_max = 1.5 ! dynamics nudging magnitude above the ABL [hour] (~1 rn_Dt)
rn_ltra_min = 4.5 ! tracers nudging magnitude inside the ABL [hour] (~3 rn_Dt)
rn_ltra_max = 1.5 ! tracers nudging magnitude above the ABL [hour] (~1 rn_Dt)
+ rn_vfac = 1. ! multiplicative factor for ocean/ice velocity
nn_amxl = 0 ! mixing length: = 0 Deardorff 80 length-scale
! = 1 length-scale based on the distance to the PBL height
! = 2 Bougeault & Lacarrere 89 length-scale
diff --git a/sette/.gitignore b/sette/.gitignore
index 1e06c5435..2bbe23174 100644
--- a/sette/.gitignore
+++ b/sette/.gitignore
@@ -1,2 +1,3 @@
+param.cfg
job_batch_template
output.sette
diff --git a/sette/param.cfg b/sette/param.default
similarity index 100%
rename from sette/param.cfg
rename to sette/param.default
diff --git a/src/ABL/ablmod.F90 b/src/ABL/ablmod.F90
index 5a4dd65ae..b5b66b8a6 100644
--- a/src/ABL/ablmod.F90
+++ b/src/ABL/ablmod.F90
@@ -82,7 +82,7 @@ CONTAINS
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du ! Cd x Du (T-point)
REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: psen ! Ch x Du
REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pevp ! Ce x Du
- REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! ||uwnd||
+ REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! ||uwnd - uoce||
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: plat ! latent heat flux
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui ! taux
REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj ! tauy
@@ -103,6 +103,8 @@ CONTAINS
#endif
!
REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j
+ REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zsspt
+ REAL(wp), DIMENSION(1:jpi,1:jpj ) :: ztabs, zpre
REAL(wp), DIMENSION(1:jpi ,2:jpka) :: zCF
!
REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_a
@@ -141,6 +143,9 @@ CONTAINS
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(:,:)+rt0, pslp_dta(:,:) )
+
!
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! 1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
@@ -186,7 +191,7 @@ CONTAINS
IF(jtra == jp_ta) THEN
DO ji = 1,jpi ! surface boundary condition for temperature
zztmp1 = psen(ji, jj)
- zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )
+ zztmp2 = psen(ji, jj) * zsspt(ji, jj)
#if defined key_si3
zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)
zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj)
@@ -340,26 +345,26 @@ CONTAINS
DO jk = 3, jpkam1
DO ji = 1, jpi
- z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
- z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
- z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
+ z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
+ z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
+ z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
END DO
END DO
- DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
+ DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
!++ Surface boundary condition
z_elem_a( ji, 2 ) = 0._wp
z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
!
- zztmp1 = pcd_du(ji, jj)
- zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) )
+ zztmp1 = pcd_du(ji, jj)
+ zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji, jj ) ) * rn_vfac
#if defined key_si3
zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
- zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj) )
+ zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj ) ) * rn_vfac
zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
#endif
z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
- u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
+ u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2
! idealized test cases only
!IF( ln_topbc_neumann ) THEN
@@ -380,11 +385,10 @@ CONTAINS
!!
!! Matrix inversion
!! ----------------------------------------------------------
- !DO ji = 2, jpi
- DO ji = 1, jpi !!GS: TBI
+ DO ji = 1, jpi !DO ji = 2, jpi !!GS: TBI
zcff = 1._wp / z_elem_b( ji, 2 )
- zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
- u_abl (ji,jj,2,nt_a) = zcff * u_abl(ji,jj,2,nt_a)
+ zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
+ u_abl (ji,jj,2,nt_a) = zcff * u_abl ( ji, jj, 2, nt_a )
END DO
DO jk = 3, jpka
@@ -414,9 +418,9 @@ CONTAINS
!
DO jk = 3, jpkam1
DO ji = 1, jpi
- z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
- z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
- z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
+ z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal
+ z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal
+ z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal
END DO
END DO
@@ -426,10 +430,10 @@ CONTAINS
z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
!
zztmp1 = pcd_du(ji, jj)
- zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) )
+ zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) * rn_vfac
#if defined key_si3
zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj)
- zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) )
+ zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) ) * rn_vfac
zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice
#endif
z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1
@@ -455,9 +459,9 @@ CONTAINS
!! Matrix inversion
!! ----------------------------------------------------------
DO ji = 1, jpi
- zcff = 1._wp / z_elem_b( ji, 2 )
- zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
- v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a )
+ zcff = 1._wp / z_elem_b( ji, 2 )
+ zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 )
+ v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a )
END DO
DO jk = 3, jpka
@@ -494,8 +498,8 @@ CONTAINS
zmsk = msk_abl(ji,jj)
zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2 &
& + jp_alp1_dyn * zsig + jp_alp0_dyn
- zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points
- ! rn_Dt = rDt_abl / nn_fsbc
+ zcff = (1._wp-zmsk) + zmsk * zcff2 * rDt_abl ! zcff = 1 for masked points
+ ! rn_Dt = rDt_abl / nn_fsbc
zcff = zcff * rest_eq(ji,jj)
u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) &
& + zcff * pu_dta( ji, jj, jk )
@@ -517,8 +521,8 @@ CONTAINS
zmsk = msk_abl(ji,jj)
zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2 &
& + jp_alp1_tra * zsig + jp_alp0_tra
- zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points
- ! rn_Dt = rDt_abl / nn_fsbc
+ zcff = (1._wp-zmsk) + zmsk * zcff2 * rDt_abl ! zcff = 1 for masked points
+ ! rn_Dt = rDt_abl / nn_fsbc
tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) &
& + zcff * pt_dta( ji, jj, jk )
@@ -534,7 +538,7 @@ CONTAINS
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
CALL lbc_lnk( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1._wp, v_abl(:,:,:,nt_a) , 'T', -1._wp )
- CALL lbc_lnk( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1._wp , tq_abl(:,:,:,nt_a,jp_qa), 'T', 1._wp , kfillmode = jpfillnothing ) ! ++++ this should not be needed...
+ !CALL lbc_lnk( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1._wp, tq_abl(:,:,:,nt_a,jp_qa), 'T', 1._wp, kfillmode = jpfillnothing ) ! ++++ this should not be needed...
!
#if defined key_xios
! 2D & first ABL level
@@ -553,8 +557,8 @@ CONTAINS
IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) )
END IF
IF( ln_hpgls_frc ) THEN
- IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
- IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) )
+ IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", -pgv_dta(:,:,2)/MAX( ABS( fft_abl(:,:)), 2.5e-5_wp)*fft_abl(:,:)/ABS(fft_abl(:,:)) )
+ IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgu_dta(:,:,2)/MAX( ABS( fft_abl(:,:)), 2.5e-5_wp)*fft_abl(:,:)/ABS(fft_abl(:,:)) )
END IF
! 3D (all ABL levels)
IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) )
@@ -576,8 +580,8 @@ CONTAINS
IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgv_dta(:,:,2:jpka) )
END IF
IF( ln_hpgls_frc ) THEN
- IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) )
- IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", -pgv_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) )
+ IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", -pgv_dta(:,:,2:jpka)/MAX( ABS( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:))), 2.5e-5_wp) * RESHAPE( fft_abl(:,:)/ABS(fft_abl(:,:)), (/jpi,jpj,jpka-1/), fft_abl(:,:)/ABS(fft_abl(:,:))) )
+ IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgu_dta(:,:,2:jpka)/MAX( ABS( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:))), 2.5e-5_wp) * RESHAPE( fft_abl(:,:)/ABS(fft_abl(:,:)), (/jpi,jpj,jpka-1/), fft_abl(:,:)/ABS(fft_abl(:,:))) )
END IF
#endif
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -585,18 +589,19 @@ CONTAINS
! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta )
- zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa )
- zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) )
- psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) !GS: negative sign to respect aerobulk convention
- pevp( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
+ ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta )
+ zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa )
+ zpre( ji, jj ) = pres_temp( zhumi, pslp_dta(ji,jj), ght_abl(2), ptpot=ztemp, pta=ztabs( ji, jj ) )
+ zcff = rho_air( ztabs( ji, jj ), zhumi, zpre( ji, jj ) )
+ psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( zsspt(ji,jj) - ztemp ) !GS: negative sign to respect aerobulk convention
+ pevp( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) )
plat( ji, jj ) = - L_vap( psst(ji,jj) ) * pevp( ji, jj )
rhoa( ji, jj ) = zcff
END_2D
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
- zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * ( pssu(ji ,jj) + pssu(ji-1,jj) )
- zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * ( pssv(ji,jj ) + pssv(ji,jj-1) )
+ zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * ( pssu(ji ,jj) + pssu(ji-1,jj) ) * rn_vfac
+ zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * ( pssv(ji,jj ) + pssv(ji,jj-1) ) * rn_vfac
END_2D
!
CALL lbc_lnk( 'ablmod', zwnd_i(:,:) , 'T', -1.0_wp, zwnd_j(:,:) , 'T', -1.0_wp )
@@ -638,16 +643,17 @@ CONTAINS
! ------------------------------------------------------------ !
! Wind stress relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
- DO_2D( 0, 0, 0, 0 )
- ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) &
- & * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
- ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) &
- & * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
- END_2D
- CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
- !
- IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' &
- & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' )
+ !DO_2D( 0, 0, 0, 0 )
+ ! ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) &
+ ! & * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
+ ! ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) &
+ ! & * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
+ !END_2D
+ !CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
+ !!
+ !IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' &
+ ! & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' )
+
! ------------------------------------------------------------ !
! Wind stress relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
@@ -658,10 +664,10 @@ CONTAINS
ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) &
& + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) &
- & * ( zztmp1 - pssu_ice(ji,jj) )
+ & * ( zztmp1 - pssu_ice(ji,jj) * rn_vfac )
ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) &
& + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) &
- & * ( zztmp2 - pssv_ice(ji,jj) )
+ & * ( zztmp2 - pssv_ice(ji,jj) * rn_vfac )
END_2D
CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice,'V', -1.0_wp )
!
diff --git a/src/ABL/par_abl.F90 b/src/ABL/par_abl.F90
index 432dda3b7..35dcda6f3 100644
--- a/src/ABL/par_abl.F90
+++ b/src/ABL/par_abl.F90
@@ -60,6 +60,7 @@ MODULE par_abl
REAL(wp), PUBLIC :: rn_ltra_min !: namelist parameter
REAL(wp), PUBLIC :: rn_ltra_max !: namelist parameter
REAL(wp), PUBLIC :: rn_Ric !: critical Richardson number
+ REAL(wp), PUBLIC :: rn_vfac !: multiplicative factor for ocean/ice velocity
!!---------------------------------------------------------------------
!! ABL parameters for the vertical profile of the restoring term
diff --git a/src/ABL/sbcabl.F90 b/src/ABL/sbcabl.F90
index 9c56d719e..f087614e3 100644
--- a/src/ABL/sbcabl.F90
+++ b/src/ABL/sbcabl.F90
@@ -71,7 +71,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, &
- & ln_smth_pblh
+ & rn_vfac, ln_smth_pblh
!!---------------------------------------------------------------------
! Namelist namsbc_abl in reference namelist : ABL parameters
@@ -217,8 +217,8 @@ CONTAINS
& + jp_alp1_dyn * jp_bmax + jp_alp0_dyn ) * rDt_abl
IF(lwp) THEN
IF(nn_dyn_restore > 0) THEN
- WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff
- WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1
+ WRITE(numout,*) ' Minimum value for ABL dynamics restoring = ',zcff
+ WRITE(numout,*) ' Maximum value for ABL dynamics restoring = ',zcff1
! Check that restoring coefficients are between 0 and 1
IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp ) &
& CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' )
@@ -235,8 +235,8 @@ CONTAINS
zcff1 = ( jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2 &
& + jp_alp1_tra * jp_bmax + jp_alp0_tra ) * rDt_abl
IF(lwp) THEN
- WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff
- WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1
+ WRITE(numout,*) ' Minimum value for ABL tracers restoring = ',zcff
+ WRITE(numout,*) ' Maximum value for ABL tracers restoring = ',zcff1
! Check that restoring coefficients are between 0 and 1
IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp ) &
& CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' )
@@ -369,6 +369,7 @@ CONTAINS
& , utau_ice, vtau_ice & ! =>> out
#endif
& )
+
!!-------------------------------------------------------------------------------------------
!! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since
!! time swap is done in abl_stp
diff --git a/src/OCE/SBC/sbc_phy.F90 b/src/OCE/SBC/sbc_phy.F90
index fd527e388..cabec70fc 100644
--- a/src/OCE/SBC/sbc_phy.F90
+++ b/src/OCE/SBC/sbc_phy.F90
@@ -22,20 +22,25 @@ MODULE sbc_phy
USE phycst ! physical constants
IMPLICIT NONE
- !PRIVATE
- PUBLIC !! Haleluja that was the solution
+ PUBLIC !! Haleluja that was the solution for AGRIF
INTEGER , PARAMETER, PUBLIC :: nb_iter0 = 5 ! Default number of itterations in bulk-param algorithms (can be overriden b.m.o `nb_iter` optional argument)
!! (mainly removed from sbcblk.F90)
- REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp !: Specic heat of dry air, constant pressure [J/K/kg]
- REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp !: Specic heat of water vapor, constant pressure [J/K/kg]
- REAL(wp), PARAMETER, PUBLIC :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
- REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
- REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
- REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
- REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...)
- REAL(wp), PARAMETER, PUBLIC :: albo = 0.066_wp !: ocean albedo assumed to be constant
+ REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp !: Specic heat of dry air, constant pressure [J/K/kg]
+ REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp !: Specic heat of water vapor, constant pressure [J/K/kg]
+ REAL(wp), PARAMETER, PUBLIC :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
+ REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
+ REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
+ REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
+ REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...)
+ REAL(wp), PARAMETER, PUBLIC :: albo = 0.066_wp !: ocean albedo assumed to be constant
+ REAL(wp), PARAMETER, PUBLIC :: R_gas = 8.314510_wp !: Universal molar gas constant [J/mol/K]
+ REAL(wp), PARAMETER, PUBLIC :: rmm_dryair = 28.9647e-3_wp !: dry air molar mass / molecular weight [kg/mol]
+ REAL(wp), PARAMETER, PUBLIC :: rmm_water = 18.0153e-3_wp !: water molar mass / molecular weight [kg/mol]
+ REAL(wp), PARAMETER, PUBLIC :: rmm_ratio = rmm_water / rmm_dryair
+ REAL(wp), PARAMETER, PUBLIC :: rgamma_dry = R_gas / ( rmm_dryair * rCp_dry ) !: Poisson constant for dry air
+ REAL(wp), PARAMETER, PUBLIC :: rpref = 1.e5_wp !: reference air pressure for exner function [Pa]
!
REAL(wp), PARAMETER, PUBLIC :: rho0_a = 1.2_wp !: Approx. of density of air [kg/m^3]
REAL(wp), PARAMETER, PUBLIC :: rho0_w = 1025._wp !: Density of sea-water (ECMWF->1025) [kg/m^3]
@@ -82,6 +87,14 @@ MODULE sbc_phy
MODULE PROCEDURE virt_temp_vctr, virt_temp_sclr
END INTERFACE virt_temp
+ INTERFACE pres_temp
+ MODULE PROCEDURE pres_temp_vctr, pres_temp_sclr
+ END INTERFACE pres_temp
+
+ INTERFACE theta_exner
+ MODULE PROCEDURE theta_exner_vctr, theta_exner_sclr
+ END INTERFACE theta_exner
+
INTERFACE visc_air
MODULE PROCEDURE visc_air_vctr, visc_air_sclr
END INTERFACE visc_air
@@ -97,6 +110,7 @@ MODULE sbc_phy
INTERFACE e_sat_ice
MODULE PROCEDURE e_sat_ice_vctr, e_sat_ice_sclr
END INTERFACE e_sat_ice
+
INTERFACE de_sat_dt_ice
MODULE PROCEDURE de_sat_dt_ice_vctr, de_sat_dt_ice_sclr
END INTERFACE de_sat_dt_ice
@@ -147,6 +161,8 @@ MODULE sbc_phy
PUBLIC virt_temp
+ PUBLIC pres_temp
+ PUBLIC theta_exner
PUBLIC rho_air
PUBLIC visc_air
PUBLIC L_vap
@@ -157,7 +173,7 @@ MODULE sbc_phy
PUBLIC q_sat
PUBLIC q_air_rh
PUBLIC dq_sat_dt_ice
- !:
+ !
PUBLIC update_qnsol_tau
PUBLIC alpha_sw
PUBLIC bulk_formula
@@ -204,14 +220,140 @@ CONTAINS
!! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
!
END FUNCTION virt_temp_sclr
- !!
+
FUNCTION virt_temp_vctr( pta, pqa )
+
REAL(wp), DIMENSION(jpi,jpj) :: virt_temp_vctr !: virtual temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air [kg/kg]
+
virt_temp_vctr(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
+
END FUNCTION virt_temp_vctr
- !===============================================================================================
+
+
+ FUNCTION pres_temp_sclr( pqspe, pslp, pz, ptpot, pta, l_ice )
+
+ !!-------------------------------------------------------------------------------
+ !! *** FUNCTION pres_temp ***
+ !!
+ !! ** Purpose : compute air pressure using barometric equation
+ !! from either potential or absolute air temperature
+ !! ** Author: G. Samson, Feb 2021
+ !!-------------------------------------------------------------------------------
+
+ REAL(wp) :: pres_temp_sclr ! air pressure [Pa]
+ REAL(wp), INTENT(in ) :: pqspe ! air specific humidity [kg/kg]
+ REAL(wp), INTENT(in ) :: pslp ! sea-level pressure [Pa]
+ REAL(wp), INTENT(in ) :: pz ! height above surface [m]
+ REAL(wp), INTENT(in ) , OPTIONAL :: ptpot ! air potential temperature [K]
+ REAL(wp), INTENT(inout), OPTIONAL :: pta ! air absolute temperature [K]
+ REAL(wp) :: ztpot, zta, zpa, zxm, zmask, zqsat
+ INTEGER :: it, niter = 3 ! iteration indice and number
+ LOGICAL , INTENT(in) , OPTIONAL :: l_ice ! sea-ice presence
+ LOGICAL :: lice ! sea-ice presence
+
+ IF( PRESENT(ptpot) ) THEN
+ zmask = 1._wp
+ ztpot = ptpot
+ zta = 0._wp
+ ELSE
+ zmask = 0._wp
+ ztpot = 0._wp
+ zta = pta
+ ENDIF
+
+ lice = .FALSE.
+ IF( PRESENT(l_ice) ) lice = l_ice
+
+ zpa = pslp ! air pressure first guess [Pa]
+ DO it = 1, niter
+ zta = ztpot * ( zpa / rpref )**rgamma_dry * zmask + (1._wp - zmask) * zta
+ zqsat = q_sat( zta, zpa, l_ice=lice ) ! saturation specific humidity [kg/kg]
+ zxm = (1._wp - pqspe/zqsat) * rmm_dryair + pqspe/zqsat * rmm_water ! moist air molar mass [kg/mol]
+ zpa = pslp * EXP( -grav * zxm * pz / ( R_gas * zta ) )
+ END DO
+
+ pres_temp_sclr = zpa
+ IF(( PRESENT(pta) ).AND.( PRESENT(ptpot) )) pta = zta
+
+ END FUNCTION pres_temp_sclr
+
+
+ FUNCTION pres_temp_vctr( pqspe, pslp, pz, ptpot, pta, l_ice )
+
+ !!-------------------------------------------------------------------------------
+ !! *** FUNCTION pres_temp ***
+ !!
+ !! ** Purpose : compute air pressure using barometric equation
+ !! from either potential or absolute air temperature
+ !! ** Author: G. Samson, Feb 2021
+ !!-------------------------------------------------------------------------------
+
+ REAL(wp), DIMENSION(jpi,jpj) :: pres_temp_vctr ! air pressure [Pa]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqspe ! air specific humidity [kg/kg]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pslp ! sea-level pressure [Pa]
+ REAL(wp), INTENT(in ) :: pz ! height above surface [m]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) , OPTIONAL :: ptpot ! air potential temperature [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(inout), OPTIONAL :: pta ! air absolute temperature [K]
+ INTEGER :: ji, jj ! loop indices
+ LOGICAL , INTENT(in) , OPTIONAL :: l_ice ! sea-ice presence
+ LOGICAL :: lice ! sea-ice presence
+
+ lice = .FALSE.
+ IF( PRESENT(l_ice) ) lice = l_ice
+
+ IF( PRESENT(ptpot) ) THEN
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, ptpot=ptpot(ji,jj), pta=pta(ji,jj), l_ice=lice )
+ END_2D
+ ELSE
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, pta=pta(ji,jj), l_ice=lice )
+ END_2D
+ ENDIF
+
+ END FUNCTION pres_temp_vctr
+
+
+ FUNCTION theta_exner_sclr( pta, ppa )
+
+ !!-------------------------------------------------------------------------------
+ !! *** FUNCTION theta_exner ***
+ !!
+ !! ** Purpose : compute air/surface potential temperature from absolute temperature
+ !! and pressure using Exner function
+ !! ** Author: G. Samson, Feb 2021
+ !!-------------------------------------------------------------------------------
+
+ REAL(wp) :: theta_exner_sclr ! air/surface potential temperature [K]
+ REAL(wp), INTENT(in) :: pta ! air/surface absolute temperature [K]
+ REAL(wp), INTENT(in) :: ppa ! air/surface pressure [Pa]
+
+ theta_exner_sclr = pta * ( rpref / ppa ) ** rgamma_dry
+
+ END FUNCTION theta_exner_sclr
+
+ FUNCTION theta_exner_vctr( pta, ppa )
+
+ !!-------------------------------------------------------------------------------
+ !! *** FUNCTION theta_exner ***
+ !!
+ !! ** Purpose : compute air/surface potential temperature from absolute temperature
+ !! and pressure using Exner function
+ !! ** Author: G. Samson, Feb 2021
+ !!-------------------------------------------------------------------------------
+
+ REAL(wp), DIMENSION(jpi,jpj) :: theta_exner_vctr ! air/surface potential temperature [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta ! air/surface absolute temperature [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! air/surface pressure [Pa]
+ INTEGER :: ji, jj ! loop indices
+
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ theta_exner_vctr(ji,jj) = theta_exner_sclr( pta(ji,jj), ppa(ji,jj) )
+ END_2D
+
+ END FUNCTION theta_exner_vctr
FUNCTION rho_air_vctr( ptak, pqa, ppa )
@@ -227,7 +369,9 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! pressure in [Pa]
REAL(wp), DIMENSION(jpi,jpj) :: rho_air_vctr ! density of moist air [kg/m^3]
!!-------------------------------------------------------------------------------
+
rho_air_vctr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
+
END FUNCTION rho_air_vctr
FUNCTION rho_air_sclr( ptak, pqa, ppa )
@@ -244,9 +388,8 @@ CONTAINS
REAL(wp) :: rho_air_sclr ! density of moist air [kg/m^3]
!!-------------------------------------------------------------------------------
rho_air_sclr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
- END FUNCTION rho_air_sclr
-
+ END FUNCTION rho_air_sclr
FUNCTION visc_air_sclr(ptak)
@@ -268,12 +411,15 @@ CONTAINS
END FUNCTION visc_air_sclr
FUNCTION visc_air_vctr(ptak)
+
REAL(wp), DIMENSION(jpi,jpj) :: visc_air_vctr ! kinetic viscosity (m^2/s)
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K)
INTEGER :: ji, jj ! dummy loop indices
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) )
+ visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) )
END_2D
+
END FUNCTION visc_air_vctr
@@ -318,10 +464,12 @@ CONTAINS
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!-------------------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
REAL(wp), DIMENSION(jpi,jpj) :: cp_air_vctr ! specific heat of moist air [J/K/kg]
!!-------------------------------------------------------------------------------
+
cp_air_vctr = rCp_dry + rCp_vap * pqa
+
END FUNCTION cp_air_vctr
FUNCTION cp_air_sclr( pqa )
@@ -335,12 +483,12 @@ CONTAINS
REAL(wp), INTENT(in) :: pqa ! air specific humidity [kg/kg]
REAL(wp) :: cp_air_sclr ! specific heat of moist air [J/K/kg]
!!-------------------------------------------------------------------------------
+
cp_air_sclr = rCp_dry + rCp_vap * pqa
- END FUNCTION cp_air_sclr
+ END FUNCTION cp_air_sclr
- !===============================================================================================
FUNCTION gamma_moist_sclr( ptak, pqa )
!!----------------------------------------------------------------------------------
!! ** Purpose : Compute the moist adiabatic lapse-rate.
@@ -365,17 +513,19 @@ CONTAINS
gamma_moist_sclr = grav * ( 1._wp + zLvap*zwa*ziRT ) / ( rCp_dry + zLvap*zLvap*zwa*reps0*ziRT/zta )
!!
END FUNCTION gamma_moist_sclr
- !!
+
FUNCTION gamma_moist_vctr( ptak, pqa )
+
REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa
INTEGER :: ji, jj
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
+ gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
END_2D
+
END FUNCTION gamma_moist_vctr
- !===============================================================================================
FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
@@ -398,17 +548,15 @@ CONTAINS
!!-------------------------------------------------------------------
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- !
- zqa = (1._wp + rctv0*pqa(ji,jj))
- !
- ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
- ! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
- ! or
- ! b/ -u* [ theta* + 0.61 theta q* ]
- !
- One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
- & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
- !
+ zqa = (1._wp + rctv0*pqa(ji,jj))
+ !
+ ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
+ ! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
+ ! or
+ ! b/ -u* [ theta* + 0.61 theta q* ]
+ !
+ One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
+ & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
END_2D
!
One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
@@ -416,7 +564,6 @@ CONTAINS
END FUNCTION One_on_L
- !===============================================================================================
FUNCTION Ri_bulk_sclr( pz, psst, ptha, pssq, pqa, pub, pta_layer, pqa_layer )
!!----------------------------------------------------------------------------------
!! Bulk Richardson number according to "wide-spread equation"...
@@ -427,7 +574,7 @@ CONTAINS
!!----------------------------------------------------------------------------------
REAL(wp) :: Ri_bulk_sclr
REAL(wp), INTENT(in) :: pz ! height above the sea (aka "delta z") [m]
- REAL(wp), INTENT(in) :: psst ! SST [K]
+ REAL(wp), INTENT(in) :: psst ! potential SST [K]
REAL(wp), INTENT(in) :: ptha ! pot. air temp. at height "pz" [K]
REAL(wp), INTENT(in) :: pssq ! 0.98*q_sat(SST) [kg/kg]
REAL(wp), INTENT(in) :: pqa ! air spec. hum. at height "pz" [kg/kg]
@@ -437,29 +584,20 @@ CONTAINS
!!
LOGICAL :: l_ptqa_l_prvd = .FALSE.
REAL(wp) :: zqa, zta, zgamma, zdthv, ztv, zsstv ! local scalars
+ REAL(wp) :: ztptv
!!-------------------------------------------------------------------
- IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE.
+ IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd = .TRUE.
!
- zsstv = virt_temp_sclr( psst, pssq ) ! virtual SST (absolute==potential because z=0!)
+ zsstv = virt_temp_sclr( psst, pssq ) ! virtual potential SST
+ ztptv = virt_temp_sclr( ptha, pqa ) ! virtual potential air temperature
+ zdthv = ztptv - zsstv ! air-sea delta of "virtual potential temperature"
!
- zdthv = virt_temp_sclr( ptha, pqa ) - zsstv ! air-sea delta of "virtual potential temperature"
- !
- !! ztv: estimate of the ABSOLUTE virtual temp. within the layer
- IF( l_ptqa_l_prvd ) THEN
- ztv = virt_temp_sclr( pta_layer, pqa_layer )
- ELSE
- zqa = 0.5_wp*( pqa + pssq ) ! ~ mean q within the layer...
- zta = 0.5_wp*( psst + ptha - gamma_moist(ptha, zqa)*pz ) ! ~ mean absolute temperature of air within the layer
- zta = 0.5_wp*( psst + ptha - gamma_moist( zta, zqa)*pz ) ! ~ mean absolute temperature of air within the layer
- zgamma = gamma_moist(zta, zqa) ! Adiabatic lapse-rate for moist air within the layer
- ztv = 0.5_wp*( zsstv + virt_temp_sclr( ptha-zgamma*pz, pqa ) )
- END IF
- !
- Ri_bulk_sclr = grav*zdthv*pz / ( ztv*pub*pub ) ! the usual definition of Ri_bulk_sclr
+ Ri_bulk_sclr = grav * zdthv * pz / ( ztptv * pub * pub ) ! the usual definition of Ri_bulk_sclr
!
END FUNCTION Ri_bulk_sclr
- !!
+
FUNCTION Ri_bulk_vctr( pz, psst, ptha, pssq, pqa, pub, pta_layer, pqa_layer )
+
REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk_vctr
REAL(wp) , INTENT(in) :: pz ! height above the sea (aka "delta z") [m]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! SST [K]
@@ -472,19 +610,22 @@ CONTAINS
!!
LOGICAL :: l_ptqa_l_prvd = .FALSE.
INTEGER :: ji, jj
- IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE.
- DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+
+ IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd = .TRUE.
IF( l_ptqa_l_prvd ) THEN
- Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), &
- & pta_layer=pta_layer(ji,jj ), pqa_layer=pqa_layer(ji,jj ) )
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), &
+ & pta_layer=pta_layer(ji,jj ), pqa_layer=pqa_layer(ji,jj ) )
+ END_2D
ELSE
- Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) )
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) )
+ END_2D
END IF
- END_2D
+
END FUNCTION Ri_bulk_vctr
- !===============================================================================================
- !===============================================================================================
+
FUNCTION e_sat_sclr( ptak )
!!----------------------------------------------------------------------------------
!! *** FUNCTION e_sat_sclr ***
@@ -509,7 +650,7 @@ CONTAINS
& + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
!
END FUNCTION e_sat_sclr
- !!
+
FUNCTION e_sat_vctr(ptak)
REAL(wp), DIMENSION(jpi,jpj) :: e_sat_vctr !: vapour pressure at saturation [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak !: temperature (K)
@@ -518,9 +659,8 @@ CONTAINS
e_sat_vctr(ji,jj) = e_sat_sclr(ptak(ji,jj))
END_2D
END FUNCTION e_sat_vctr
- !===============================================================================================
- !===============================================================================================
+
FUNCTION e_sat_ice_sclr(ptak)
!!---------------------------------------------------------------------------------
!! Same as "e_sat" but over ice rather than water!
@@ -536,8 +676,9 @@ CONTAINS
zle = rAg_i*(ztmp - 1._wp) + rBg_i*LOG10(ztmp) + rCg_i*(1._wp - zta/rtt0) + rDg_i
!!
e_sat_ice_sclr = 100._wp * 10._wp**zle
+
END FUNCTION e_sat_ice_sclr
- !!
+
FUNCTION e_sat_ice_vctr(ptak)
!! Same as "e_sat" but over ice rather than water!
REAL(wp), DIMENSION(jpi,jpj) :: e_sat_ice_vctr !: vapour pressure at saturation in presence of ice [Pa]
@@ -545,10 +686,12 @@ CONTAINS
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) )
+ e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) )
END_2D
+
END FUNCTION e_sat_ice_vctr
- !!
+
+
FUNCTION de_sat_dt_ice_sclr(ptak)
!!---------------------------------------------------------------------------------
!! d [ e_sat_ice ] / dT (derivative / temperature)
@@ -566,7 +709,7 @@ CONTAINS
!!
de_sat_dt_ice_sclr = LOG(10._wp) * zde * e_sat_ice_sclr(zta)
END FUNCTION de_sat_dt_ice_sclr
- !!
+
FUNCTION de_sat_dt_ice_vctr(ptak)
!! Same as "e_sat" but over ice rather than water!
REAL(wp), DIMENSION(jpi,jpj) :: de_sat_dt_ice_vctr !: vapour pressure at saturation in presence of ice [Pa]
@@ -574,15 +717,12 @@ CONTAINS
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) )
+ de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) )
END_2D
- END FUNCTION de_sat_dt_ice_vctr
-
+ END FUNCTION de_sat_dt_ice_vctr
- !===============================================================================================
- !===============================================================================================
FUNCTION q_sat_sclr( pta, ppa, l_ice )
!!---------------------------------------------------------------------------------
!! *** FUNCTION q_sat_sclr ***
@@ -606,9 +746,11 @@ CONTAINS
ze_s = e_sat( pta ) ! Vapour pressure at saturation (Goff) :
END IF
q_sat_sclr = reps0*ze_s/(ppa - (1._wp - reps0)*ze_s)
+
END FUNCTION q_sat_sclr
- !!
+
FUNCTION q_sat_vctr( pta, ppa, l_ice )
+
REAL(wp), DIMENSION(jpi,jpj) :: q_sat_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute temperature of air [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa !: atmospheric pressure [Pa]
@@ -619,13 +761,12 @@ CONTAINS
lice = .FALSE.
IF( PRESENT(l_ice) ) lice = l_ice
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice )
+ q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice )
END_2D
+
END FUNCTION q_sat_vctr
- !===============================================================================================
- !===============================================================================================
FUNCTION dq_sat_dt_ice_sclr( pta, ppa )
!!---------------------------------------------------------------------------------
!! *** FUNCTION dq_sat_dt_ice_sclr ***
@@ -646,21 +787,21 @@ CONTAINS
dq_sat_dt_ice_sclr = reps0*ppa*zde_s_dt / ( ztmp*ztmp )
!
END FUNCTION dq_sat_dt_ice_sclr
- !!
+
FUNCTION dq_sat_dt_ice_vctr( pta, ppa )
+
REAL(wp), DIMENSION(jpi,jpj) :: dq_sat_dt_ice_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute temperature of air [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa !: atmospheric pressure [Pa]
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) )
+ dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) )
END_2D
+
END FUNCTION dq_sat_dt_ice_vctr
- !===============================================================================================
- !===============================================================================================
FUNCTION q_air_rh(prha, ptak, ppa)
!!----------------------------------------------------------------------------------
!! Specific humidity of air out of Relative Humidity
@@ -677,14 +818,14 @@ CONTAINS
!!----------------------------------------------------------------------------------
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
- q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze)
+ ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
+ q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze)
END_2D
!
END FUNCTION q_air_rh
- SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, &
+ SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, prhoa, &
& pQns, pTau, &
& Qlat)
!!----------------------------------------------------------------------------------
@@ -704,6 +845,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prlw ! downwelling longwave radiative flux [W/m^2]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prhoa ! air density [kg/m3]
!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
@@ -714,50 +856,53 @@ CONTAINS
INTEGER :: ji, jj ! dummy loop indices
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- zdt = pTa(ji,jj) - pTs(ji,jj) ; zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
- zdq = pqa(ji,jj) - pqs(ji,jj) ; zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
- zz0 = pust(ji,jj)/pUb(ji,jj)
- zCd = zz0*zz0
- zCh = zz0*ptst(ji,jj)/zdt
- zCe = zz0*pqst(ji,jj)/zdq
- CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
- & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), &
- & pTau(ji,jj), zQsen, zQlat )
+ zdt = pTa(ji,jj) - pTs(ji,jj) ; zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
+ zdq = pqa(ji,jj) - pqs(ji,jj) ; zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
+ zz0 = pust(ji,jj)/pUb(ji,jj)
+ zCd = zz0*zz0
+ zCh = zz0*ptst(ji,jj)/zdt
+ zCe = zz0*pqst(ji,jj)/zdq
- zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux
+ CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
+ & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), &
+ & pTau(ji,jj), zQsen, zQlat )
- pQns(ji,jj) = zQlat + zQsen + zQlw
+ zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux
+
+ pQns(ji,jj) = zQlat + zQsen + zQlw
+
+ IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
- IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
END_2D
+
END SUBROUTINE UPDATE_QNSOL_TAU
SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
& pCd, pCh, pCe, &
- & pwnd, pUb, ppa, &
+ & pwnd, pUb, ppa, prhoa, &
& pTau, pQsen, pQlat, &
- & pEvap, prhoa, pfact_evap )
+ & pEvap, pfact_evap )
!!----------------------------------------------------------------------------------
- REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
- REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
- REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
- REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
- REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
+ REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
+ REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
+ REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
+ REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
+ REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
REAL(wp), INTENT(in) :: pCd
REAL(wp), INTENT(in) :: pCh
REAL(wp), INTENT(in) :: pCe
- REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
- REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
- REAL(wp), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
+ REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
+ REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
+ REAL(wp), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
+ REAL(wp), INTENT(in) :: prhoa ! Air density at z=pzu [kg/m^3]
!!
REAL(wp), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
REAL(wp), INTENT(out) :: pQsen ! [W/m^2]
REAL(wp), INTENT(out) :: pQlat ! [W/m^2]
!!
REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
- REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
!!
REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
@@ -766,16 +911,7 @@ CONTAINS
zfact_evap = 1._wp
IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
- !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
- ztaa = pTa ! first guess...
- DO jq = 1, 4
- zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa ) !#LB: why not "0.5*(pqs+pqa)" rather then "pqa" ???
- ztaa = pTa - zgamma*pzu ! Absolute temp. is slightly colder...
- END DO
- zrho = rho_air(ztaa, pqa, ppa)
- zrho = rho_air(ztaa, pqa, ppa-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
-
- zUrho = pUb*MAX(zrho, 1._wp) ! rho*U10
+ zUrho = pUb*MAX(prhoa, 1._wp) ! rho*U10
pTau = zUrho * pCd * pwnd ! Wind stress module
@@ -784,34 +920,33 @@ CONTAINS
pQlat = L_vap(pTs) * zevap
IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap
- IF( PRESENT(prhoa) ) prhoa = zrho
END SUBROUTINE BULK_FORMULA_SCLR
- !!
+
SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, &
& pCd, pCh, pCe, &
- & pwnd, pUb, ppa, &
+ & pwnd, pUb, ppa, prhoa, &
& pTau, pQsen, pQlat, &
- & pEvap, prhoa, pfact_evap )
+ & pEvap, pfact_evap )
!!----------------------------------------------------------------------------------
- REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
+ REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCh
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCe
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
+ REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prhoa ! Air density at z=pzu [kg/m^3]
!!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen ! [W/m^2]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat ! [W/m^2]
!!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
- REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
!!
REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
@@ -822,15 +957,15 @@ CONTAINS
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
- & pCd(ji,jj), pCh(ji,jj), pCe(ji,jj), &
- & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), &
- & pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj), &
- & pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap )
+ CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
+ & pCd(ji,jj), pCh(ji,jj), pCe(ji,jj), &
+ & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), &
+ & pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj), &
+ & pEvap=zevap, pfact_evap=zfact_evap )
- IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
- IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
+ IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
END_2D
+
END SUBROUTINE BULK_FORMULA_VCTR
@@ -846,6 +981,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]
!!----------------------------------------------------------------------------------
alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
+
END FUNCTION alpha_sw_vctr
FUNCTION alpha_sw_sclr( psst )
@@ -860,10 +996,10 @@ CONTAINS
REAL(wp), INTENT(in) :: psst ! sea-water temperature [K]
!!----------------------------------------------------------------------------------
alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
+
END FUNCTION alpha_sw_sclr
- !===============================================================================================
FUNCTION qlw_net_sclr( pdwlw, pts, l_ice )
!!---------------------------------------------------------------------------------
!! *** FUNCTION qlw_net_sclr ***
@@ -886,9 +1022,11 @@ CONTAINS
END IF
zt2 = pts*pts
qlw_net_sclr = zemiss*( pdwlw - stefan*zt2*zt2) ! zemiss used both as the IR albedo and IR emissivity...
+
END FUNCTION qlw_net_sclr
- !!
+
FUNCTION qlw_net_vctr( pdwlw, pts, l_ice )
+
REAL(wp), DIMENSION(jpi,jpj) :: qlw_net_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdwlw !: downwelling longwave (aka infrared, aka thermal) radiation [W/m^2]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts !: surface temperature [K]
@@ -899,12 +1037,14 @@ CONTAINS
lice = .FALSE.
IF( PRESENT(l_ice) ) lice = l_ice
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice )
+ qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice )
END_2D
+
END FUNCTION qlw_net_vctr
- !===============================================================================================
+
FUNCTION z0_from_Cd( pzu, pCd, ppsi )
+
REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd !: roughness length [m]
REAL(wp) , INTENT(in) :: pzu !: reference height zu [m]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd !: (neutral or non-neutral) drag coefficient []
@@ -920,9 +1060,12 @@ CONTAINS
!! Cd provided is the neutral-stability Cd, aka CdN :
z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) ) !LB: ok, double-checked!
END IF
+
END FUNCTION z0_from_Cd
+
FUNCTION Cd_from_z0( pzu, pz0, ppsi )
+
REAL(wp), DIMENSION(jpi,jpj) :: Cd_from_z0 !: (neutral or non-neutral) drag coefficient []
REAL(wp) , INTENT(in) :: pzu !: reference height zu [m]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0 !: roughness length [m]
@@ -939,6 +1082,7 @@ CONTAINS
Cd_from_z0 = 1._wp / LOG( pzu / pz0(:,:) )
END IF
Cd_from_z0 = vkarmn2 * Cd_from_z0 * Cd_from_z0
+
END FUNCTION Cd_from_z0
@@ -964,17 +1108,20 @@ CONTAINS
& + zstab * 1._wp / ( 1._wp + ram_louis * zts ) ! Stable Eq.(A7)
!
END FUNCTION f_m_louis_sclr
- !!
+
FUNCTION f_m_louis_vctr( pzu, pRib, pCdn, pz0 )
+
REAL(wp), DIMENSION(jpi,jpj) :: f_m_louis_vctr
REAL(wp), INTENT(in) :: pzu
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pRib
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCdn
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0
INTEGER :: ji, jj
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) )
+ f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) )
END_2D
+
END FUNCTION f_m_louis_vctr
@@ -1000,19 +1147,23 @@ CONTAINS
& + zstab * 1._wp / ( 1._wp + rah_louis * zts ) ! Stable Eq.(A7) !#LB: in paper it's "ram_louis" and not "rah_louis" typo or what????
!
END FUNCTION f_h_louis_sclr
- !!
+
FUNCTION f_h_louis_vctr( pzu, pRib, pChn, pz0 )
+
REAL(wp), DIMENSION(jpi,jpj) :: f_h_louis_vctr
REAL(wp), INTENT(in) :: pzu
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pRib
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pChn
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0
INTEGER :: ji, jj
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) )
+ f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) )
END_2D
+
END FUNCTION f_h_louis_vctr
+
FUNCTION UN10_from_ustar( pzu, pUzu, pus, ppsi )
!!----------------------------------------------------------------------------------
!! Provides the neutral-stability wind speed at 10 m
@@ -1122,5 +1273,5 @@ CONTAINS
END FUNCTION z0tq_LKB
- !!======================================================================
+
END MODULE sbc_phy
diff --git a/src/OCE/SBC/sbcblk.F90 b/src/OCE/SBC/sbcblk.F90
index 677005d82..a88f532dc 100644
--- a/src/OCE/SBC/sbcblk.F90
+++ b/src/OCE/SBC/sbcblk.F90
@@ -79,7 +79,7 @@ MODULE sbcblk
INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point
INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point
INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin)
- INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % )
+ INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity (kg/kg)
INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2)
INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2)
INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s)
@@ -170,7 +170,7 @@ MODULE sbcblk
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
- !! $Id: sbcblk.F90 14718 2021-04-16 09:43:50Z clem $
+ !! $Id: sbcblk.F90 15551 2021-11-28 20:19:36Z gsamson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
@@ -274,8 +274,8 @@ CONTAINS
& CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' )
IF( ln_ANDREAS ) &
& CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with ANDREAS algorithm' )
- IF( nn_fsbc /= 1 ) &
- & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
+ !IF( nn_fsbc /= 1 ) &
+ ! & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
END IF
IF( ln_skin_wl ) THEN
@@ -502,7 +502,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step
!!----------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp
+ REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp, zpre, ztheta
REAL(wp) :: ztst
LOGICAL :: llerr
!!----------------------------------------------------------------------
@@ -539,7 +539,7 @@ CONTAINS
! ! compute the surface ocean fluxes using bulk formulea
IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
- ! Specific humidity of air at z=rn_zqt !
+ ! Specific humidity of air at z=rn_zqt
SELECT CASE( nhumi )
CASE( np_humi_sph )
q_air_zt(:,:) = sf(jp_humi )%fnow(:,:,1) ! what we read in file is already a spec. humidity!
@@ -552,16 +552,15 @@ CONTAINS
& sf(jp_tair )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) !#LB: 0.01 => RH is % percent in file
END SELECT
- ! POTENTIAL temperature of air at z=rn_zqt
- ! based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
- ! (most reanalysis products provide absolute temp., not potential temp.)
+ ! Potential temperature of air at z=rn_zqt (most reanalysis products provide absolute temp., not potential temp.)
IF( ln_tair_pot ) THEN
! temperature read into file is already potential temperature, do nothing...
theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)
ELSE
! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...)
IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!'
- theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt
+ zpre(:,:) = pres_temp( q_air_zt(:,:), sf(jp_slp)%fnow(:,:,1), rn_zu, pta=sf(jp_tair)%fnow(:,:,1) )
+ theta_air_zt(:,:) = theta_exner( sf(jp_tair)%fnow(:,:,1), zpre(:,:) )
ENDIF
!
CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1), & ! <<= in
@@ -585,15 +584,12 @@ CONTAINS
ELSE
qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1)
ENDIF
- tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature instead ????
- !tatm_ice(:,:) = theta_air_zt(:,:) !#LB: THIS! ?
-
- qatm_ice(:,:) = q_air_zt(:,:) !#LB:
-
- tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
- sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
- wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
- wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
+ tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature (theta_air_zt) instead ????
+ qatm_ice(:,:) = q_air_zt(:,:)
+ tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac
+ sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac
+ wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1)
+ wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1)
ENDIF
#endif
!
@@ -651,6 +647,8 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean
REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean
REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean
+ REAL(wp), DIMENSION(jpi,jpj) :: zsspt ! potential sea-surface temperature [K]
+ REAL(wp), DIMENSION(jpi,jpj) :: zpre, ztabs ! air pressure [Pa] & absolute temperature [K]
REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2
!!---------------------------------------------------------------------
!
@@ -658,6 +656,9 @@ CONTAINS
! ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K)
ptsk(:,:) = pst(:,:) + rt0 ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!)
+ ! sea surface potential temperature [K]
+ zsspt(:,:) = theta_exner( ptsk(:,:), pslp(:,:) )
+
! --- cloud cover --- !
cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
@@ -704,7 +705,7 @@ CONTAINS
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
!! Backup "bulk SST" and associated spec. hum.
- zztmp1(:,:) = ptsk(:,:)
+ zztmp1(:,:) = zsspt(:,:)
zztmp2(:,:) = pssq(:,:)
ENDIF
@@ -713,33 +714,33 @@ CONTAINS
SELECT CASE( nblk )
CASE( np_NCAR )
- CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_ncar ( rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& nb_iter=nn_iter_algo )
!
CASE( np_COARE_3p0 )
- CALL turb_coare3p0( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_coare3p0( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_COARE_3p6 )
- CALL turb_coare3p6( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_coare3p6( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_ECMWF )
- CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& ln_skin_cs, ln_skin_wl, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, &
& nb_iter=nn_iter_algo, &
& Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
!
CASE( np_ANDREAS )
- CALL turb_andreas ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, &
+ CALL turb_andreas ( rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
& zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
& nb_iter=nn_iter_algo )
!
@@ -760,34 +761,43 @@ CONTAINS
IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
- !! ptsk and pssq have been updated!!!
- !!
- !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq:
+ !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zsspt, pssq & ptsk:
WHERE ( fr_i(:,:) > 0.001_wp )
! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
- ptsk(:,:) = zztmp1(:,:)
- pssq(:,:) = zztmp2(:,:)
+ zsspt(:,:) = zztmp1(:,:)
+ pssq(:,:) = zztmp2(:,:)
END WHERE
+ ! apply potential temperature increment to abolute SST
+ ptsk(:,:) = ptsk(:,:) + ( zsspt(:,:) - zztmp1(:,:) )
END IF
! Turbulent fluxes over ocean => BULK_FORMULA @ sbc_phy.F90
! -------------------------------------------------------------
IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zztmp = zU_zu(ji,jj)
wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod
pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
psen(ji,jj) = zztmp * zch_oce(ji,jj)
pevp(ji,jj) = zztmp * zce_oce(ji,jj)
- rhoa(ji,jj) = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) )
+ zpre(ji,jj) = pres_temp( pqair(ji,jj), pslp(ji,jj), rn_zu, ptpot=ptair(ji,jj), pta=ztabs(ji,jj) )
+ rhoa(ji,jj) = rho_air( ztabs(ji,jj), pqair(ji,jj), zpre(ji,jj) )
END_2D
+
ELSE !== BLK formulation ==! turbulent fluxes computation
- CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
+
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ zpre(ji,jj) = pres_temp( q_zu(ji,jj), pslp(ji,jj), rn_zu, ptpot=theta_zu(ji,jj), pta=ztabs(ji,jj) )
+ rhoa(ji,jj) = rho_air( ztabs(ji,jj), q_zu(ji,jj), zpre(ji,jj) )
+ END_2D
+
+ CALL BULK_FORMULA( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
& zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &
- & wndm(:,:), zU_zu(:,:), pslp(:,:), &
+ & wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:), &
& taum(:,:), psen(:,:), plat(:,:), &
- & pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac )
+ & pEvap=pevp(:,:), pfact_evap=rn_efac )
psen(:,:) = psen(:,:) * tmask(:,:,1)
plat(:,:) = plat(:,:) * tmask(:,:,1)
@@ -795,19 +805,19 @@ CONTAINS
pevp(:,:) = pevp(:,:) * tmask(:,:,1)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- IF( wndm(ji,jj) > 0._wp ) THEN
- zztmp = taum(ji,jj) / wndm(ji,jj)
+ IF( wndm(ji,jj) > 0._wp ) THEN
+ zztmp = taum(ji,jj) / wndm(ji,jj)
#if defined key_cyclone
- ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
- ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+ 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)
+ ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
+ ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
#endif
- ELSE
- ztau_i(ji,jj) = 0._wp
- ztau_j(ji,jj) = 0._wp
- ENDIF
+ ELSE
+ ztau_i(ji,jj) = 0._wp
+ ztau_j(ji,jj) = 0._wp
+ ENDIF
END_2D
IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
@@ -850,13 +860,13 @@ CONTAINS
CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ')
ENDIF
!
- ENDIF !IF( ln_abl )
+ ENDIF ! ln_blk / ln_abl
ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius
IF( ln_skin_cs .OR. ln_skin_wl ) THEN
CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius
- CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference...
+ CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference
ENDIF
!
END SUBROUTINE blk_oce_1
@@ -975,8 +985,8 @@ CONTAINS
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s]
- REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s]
- REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric wind at T-point [m/s]
+ REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric potential temperature at T-point [K]
+ REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric specific humidity at T-point [kg/kg]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s]
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! "
REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K]
@@ -990,11 +1000,9 @@ CONTAINS
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zootm_su ! sea-ice surface mean temperature
REAL(wp) :: zztmp1, zztmp2 ! temporary scalars
- REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! temporary array
+ REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
!!---------------------------------------------------------------------
!
- ! LB: ptsui is in K !!!
- !
! ------------------------------------------------------------ !
! Wind module relative to the moving ice ( U10m - U_ice ) !
! ------------------------------------------------------------ !
@@ -1003,9 +1011,10 @@ CONTAINS
wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
END_2D
!
- ! Make ice-atm. drag dependent on ice concentration
-
+ ! potential sea-ice surface temperature [K]
+ zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
+ ! sea-ice <-> atmosphere bulk transfer coefficients
SELECT CASE( nblk_ice )
CASE( np_ice_cst )
@@ -1019,17 +1028,17 @@ CONTAINS
CASE( np_ice_an05 ) ! calculate new drag from Lupkes(2015) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
- CALL turb_ice_an05( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, &
+ CALL turb_ice_an05( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
CASE( np_ice_lu12 )
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
- CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
+ CALL turb_ice_lu12( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations
ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
- CALL turb_ice_lg15( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, &
+ CALL turb_ice_lg15( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
& Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
!!
END SELECT
@@ -1049,8 +1058,8 @@ CONTAINS
! supress moving ice in wind stress computation as we don't know how to do it properly...
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)
+ pvtaui(ji,jj) = zztmp1 * pwndj(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 !!!
@@ -1073,15 +1082,15 @@ CONTAINS
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' &
& , tab2d_2=pvtaui , clinfo2=' pvtaui : ' )
ELSE ! ln_abl
+
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj)
- pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
- pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
+ pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj)
+ pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
+ pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
END_2D
- !#LB:
pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq !
- !#LB.
- ENDIF !IF( ln_blk )
+
+ ENDIF ! ln_blk / ln_abl
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ')
!
@@ -1112,7 +1121,7 @@ CONTAINS
REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow
!!
INTEGER :: ji, jj, jl ! dummy loop indices
- REAL(wp) :: zst, zst3, zsq ! local variable
+ REAL(wp) :: zst, zst3, zsq, zsipt ! local variable
REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - -
REAL(wp) :: zztmp, zzblk, zztmp1, z1_rLsub ! - -
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice
@@ -1120,13 +1129,12 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice
REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3)
+ REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2
REAL(wp), DIMENSION(jpi,jpj) :: ztri
REAL(wp), DIMENSION(jpi,jpj) :: zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg)
!!---------------------------------------------------------------------
!
zcoef_dqlw = 4._wp * emiss_i * stefan ! local scalars
- !
-
zztmp = 1. / ( 1. - albo )
dqla_ice(:,:,:) = 0._wp
@@ -1140,52 +1148,53 @@ CONTAINS
! ! ========================== !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
- zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K]
- zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present
-
- ! ----------------------------!
- ! I Radiative FLUXES !
- ! ----------------------------!
- ! Short Wave (sw)
- qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
-
- ! Long Wave (lw)
- zst3 = zst * zst * zst
- z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1)
- ! lw sensitivity
- z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
-
- ! ----------------------------!
- ! II Turbulent FLUXES !
- ! ----------------------------!
-
- ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
-
- ! Common term in bulk F. equations...
- zzblk = rhoa(ji,jj) * wndm_ice(ji,jj)
-
- ! Sensible Heat
- zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj)
- z_qsb (ji,jj,jl) = zztmp1 * (zst - theta_zu_i(ji,jj))
- z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice)
-
- ! Latent Heat
- zztmp1 = zzblk * rLsub * Ce_ice(ji,jj)
- qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ???
- IF(qla_ice(ji,jj,jl)>0._wp) dqla_ice(ji,jj,jl) = zztmp1*dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
- ! !#LB: dq_sat_dt_ice() in "sbc_phy.F90"
- !#LB: without this unjustified "condensation sensure":
- !qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj))
- !dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
-
-
- ! ----------------------------!
- ! III Total FLUXES !
- ! ----------------------------!
- ! Downward Non Solar flux
- qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
- ! Total non solar heat flux sensitivity for ice
- dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ????
+ zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K]
+ zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present
+ zsipt = theta_exner( zst, pslp(ji,jj) ) ! potential sea-ice surface temperature [K]
+
+ ! ----------------------------!
+ ! I Radiative FLUXES !
+ ! ----------------------------!
+ ! Short Wave (sw)
+ qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
+
+ ! Long Wave (lw)
+ zst3 = zst * zst * zst
+ z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1)
+ ! lw sensitivity
+ z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
+
+ ! ----------------------------!
+ ! II Turbulent FLUXES !
+ ! ----------------------------!
+
+ ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
+
+ ! Common term in bulk F. equations...
+ zzblk = rhoa(ji,jj) * wndm_ice(ji,jj)
+
+ ! Sensible Heat
+ zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj)
+ z_qsb (ji,jj,jl) = zztmp1 * (zsipt - theta_zu_i(ji,jj))
+ z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice)
+
+ ! Latent Heat
+ zztmp1 = zzblk * rLsub * Ce_ice(ji,jj)
+ qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ???
+ IF(qla_ice(ji,jj,jl)>0._wp) dqla_ice(ji,jj,jl) = zztmp1*dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
+ ! !#LB: dq_sat_dt_ice() in "sbc_phy.F90"
+ !#LB: without this unjustified "condensation sensure":
+ !qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj))
+ !dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT)
+
+
+ ! ----------------------------!
+ ! III Total FLUXES !
+ ! ----------------------------!
+ ! Downward Non Solar flux
+ qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
+ ! Total non solar heat flux sensitivity for ice
+ dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ????
END_2D
!
diff --git a/src/OCE/SBC/sbcblk_algo_coare3p0.F90 b/src/OCE/SBC/sbcblk_algo_coare3p0.F90
index 1dfd38e8d..2ad244fcb 100644
--- a/src/OCE/SBC/sbcblk_algo_coare3p0.F90
+++ b/src/OCE/SBC/sbcblk_algo_coare3p0.F90
@@ -188,6 +188,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t
REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu
REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
+ REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
!
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t ! stability parameter at height zt
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST
@@ -320,11 +321,16 @@ CONTAINS
q_zu = q_zt - q_star/vkarmn*ztmp1
ENDIF
+ IF(( l_use_cs ).OR.( l_use_wl )) THEN
+ zpre(:,:) = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) )
+ zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) )
+ ENDIF
+
IF( l_use_cs ) THEN
!! Cool-skin contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
- & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
+ & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u
CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 ) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
@@ -335,7 +341,7 @@ CONTAINS
IF( l_use_wl ) THEN
!! Warm-layer contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
& ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u
!! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) )
diff --git a/src/OCE/SBC/sbcblk_algo_coare3p6.F90 b/src/OCE/SBC/sbcblk_algo_coare3p6.F90
index 8255e81c0..cb7fff12b 100644
--- a/src/OCE/SBC/sbcblk_algo_coare3p6.F90
+++ b/src/OCE/SBC/sbcblk_algo_coare3p6.F90
@@ -178,6 +178,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t
REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu
REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
+ REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
!
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t ! stability parameter at height zt
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST
@@ -310,11 +311,16 @@ CONTAINS
q_zu = q_zt - q_star/vkarmn*ztmp1
ENDIF
+ IF(( l_use_cs ).OR.( l_use_wl )) THEN
+ zpre(:,:) = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) )
+ zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) )
+ ENDIF
+
IF( l_use_cs ) THEN
!! Cool-skin contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
- & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
+ & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u
CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 ) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
@@ -325,7 +331,7 @@ CONTAINS
IF( l_use_wl ) THEN
!! Warm-layer contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
& ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u
!! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) )
diff --git a/src/OCE/SBC/sbcblk_algo_ecmwf.F90 b/src/OCE/SBC/sbcblk_algo_ecmwf.F90
index 7d1c88631..d86c99756 100644
--- a/src/OCE/SBC/sbcblk_algo_ecmwf.F90
+++ b/src/OCE/SBC/sbcblk_algo_ecmwf.F90
@@ -183,6 +183,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air
REAL(wp), DIMENSION(jpi,jpj) :: Linv !: 1/L (inverse of Monin Obukhov length...
REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q
+ REAL(wp), DIMENSION(jpi,jpj) :: zrhoa, zpre, zta ! air pressure [Pa], density [kg/m3] & absolute temperature [k]
!
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST
!
@@ -346,12 +347,16 @@ CONTAINS
func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
+ IF(( l_use_cs ).OR.( l_use_wl )) THEN
+ zpre(:,:) = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) )
+ zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) )
+ ENDIF
IF( l_use_cs ) THEN
!! Cool-skin contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
- & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
+ & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0
CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst ) ! Qnsol -> ztmp1
@@ -363,7 +368,7 @@ CONTAINS
IF( l_use_wl ) THEN
!! Warm-layer contribution
- CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
+ CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, &
& ztmp1, ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp2
CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst )
!! Updating T_s and q_s !!!
diff --git a/src/OCE/SBC/sbccpl.F90 b/src/OCE/SBC/sbccpl.F90
index d6e020564..05070b410 100644
--- a/src/OCE/SBC/sbccpl.F90
+++ b/src/OCE/SBC/sbccpl.F90
@@ -53,7 +53,7 @@ MODULE sbccpl
USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut
#endif
- USE sbc_phy, ONLY : pp_cldf
+ USE sbc_phy, ONLY : pp_cldf, rpref
IMPLICIT NONE
PRIVATE
@@ -219,9 +219,6 @@ MODULE sbccpl
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_last_couple !: Ice fractional area at last coupling time
#endif
- REAL(wp) :: rpref = 101000._wp ! reference atmospheric pressure[N/m2]
- REAL(wp) :: r1_grau ! = 1.e0 / (grav * rho0)
-
INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) :: nrcvinfo ! OASIS info argument
!! Substitution
@@ -229,7 +226,7 @@ MODULE sbccpl
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
- !! $Id: sbccpl.F90 15004 2021-06-16 10:33:18Z mathiot $
+ !! $Id: sbccpl.F90 15551 2021-11-28 20:19:36Z gsamson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
@@ -1188,6 +1185,7 @@ CONTAINS
REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
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
!!----------------------------------------------------------------------
!
--
GitLab