Changeset 13214 for NEMO/trunk/src/ABL
- Timestamp:
- 2020-07-02T11:09:01+02:00 (4 years ago)
- Location:
- NEMO/trunk/src/ABL
- Files:
-
- 5 edited
-
abl.F90 (modified) (2 diffs)
-
ablmod.F90 (modified) (31 diffs)
-
ablrst.F90 (modified) (3 diffs)
-
par_abl.F90 (modified) (2 diffs)
-
sbcabl.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ABL/abl.F90
r12489 r13214 29 29 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: avm_abl !: turbulent viscosity [m2/s] 30 30 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: avt_abl !: turbulent diffusivity [m2/s] 31 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: mxl_abl !: mixing length [m] 31 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: mxld_abl !: dissipative mixing length [m] 32 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: mxlm_abl !: master mixing length [m] 32 33 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: tke_abl !: turbulent kinetic energy [m2/s2] 33 34 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: fft_abl !: Coriolis parameter [1/s] … … 55 56 !!---------------------------------------------------------------------- 56 57 ! 57 ALLOCATE( u_abl (1:jpi,1:jpj,1:jpka,jptime), & 58 & v_abl (1:jpi,1:jpj,1:jpka,jptime), & 59 & tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), & 60 & avm_abl(1:jpi,1:jpj,1:jpka), & 61 & avt_abl(1:jpi,1:jpj,1:jpka), & 62 & mxl_abl(1:jpi,1:jpj,1:jpka), & 63 & tke_abl(1:jpi,1:jpj,1:jpka,jptime), & 64 & fft_abl(1:jpi,1:jpj), & 65 & pblh (1:jpi,1:jpj), & 66 & msk_abl(1:jpi,1:jpj), & 67 & rest_eq(1:jpi,1:jpj), & 68 & e3t_abl(1:jpka), e3w_abl(1:jpka), ght_abl(1:jpka), ghw_abl(1:jpka), STAT=ierr ) 58 ALLOCATE( u_abl (1:jpi,1:jpj,1:jpka,jptime ), & 59 & v_abl (1:jpi,1:jpj,1:jpka,jptime ), & 60 & tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), & 61 & tke_abl (1:jpi,1:jpj,1:jpka,jptime ), & 62 & avm_abl (1:jpi,1:jpj,1:jpka ), & 63 & avt_abl (1:jpi,1:jpj,1:jpka ), & 64 & mxld_abl(1:jpi,1:jpj,1:jpka ), & 65 & mxlm_abl(1:jpi,1:jpj,1:jpka ), & 66 & fft_abl (1:jpi,1:jpj ), & 67 & pblh (1:jpi,1:jpj ), & 68 & msk_abl (1:jpi,1:jpj ), & 69 & rest_eq (1:jpi,1:jpj ), & 70 & e3t_abl (1:jpka), e3w_abl(1:jpka) , & 71 & ght_abl (1:jpka), ghw_abl(1:jpka) , STAT=ierr ) 69 72 ! 70 73 abl_alloc = ierr -
NEMO/trunk/src/ABL/ablmod.F90
r13208 r13214 2 2 !!====================================================================== 3 3 !! *** MODULE ablmod *** 4 !! Surface module : ABL computation to provide atmospheric data 4 !! Surface module : ABL computation to provide atmospheric data 5 5 !! for surface fluxes computation 6 6 !!====================================================================== 7 7 !! History : 3.6 ! 2019-03 (F. Lemarié & G. Samson) Original code 8 8 !!---------------------------------------------------------------------- 9 9 10 10 !!---------------------------------------------------------------------- 11 11 !! abl_stp : ABL single column model … … 16 16 17 17 USE phycst ! physical constants 18 USE dom_oce, ONLY : tmask 18 USE dom_oce, ONLY : tmask 19 19 USE sbc_oce, ONLY : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa 20 USE sbcblk , ONLY : rn_efac20 USE sbcblk ! use rn_efac, cdn_oce 21 21 USE sbcblk_phy ! use some physical constants for flux computation 22 22 ! … … 30 30 31 31 PUBLIC abl_stp ! called by sbcabl.F90 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustar2, zrough 33 33 !! * Substitutions 34 34 # include "do_loop_substitute.h90" … … 38 38 39 39 !=================================================================================================== 40 SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &! in41 & pu_dta, pv_dta, pt_dta, pq_dta, & 40 SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, & ! in 41 & pu_dta, pv_dta, pt_dta, pq_dta, & 42 42 & pslp_dta, pgu_dta, pgv_dta, & 43 & pcd_du, psen, pevp, & ! in/out 44 & pwndm, ptaui, ptauj, ptaum & 45 #if defined key_si3 46 & , ptm_su,pssu_ice,pssv_ice,pssq_ice,pcd_du_ice & 47 & , psen_ice, pevp_ice, pwndm_ice, pfrac_oce & 48 & , ptaui_ice, ptauj_ice & 49 #endif 50 & ) 43 & pcd_du, psen, pevp, & ! in/out 44 & pwndm, ptaui, ptauj, ptaum & 45 #if defined key_si3 46 & , ptm_su, pssu_ice, pssv_ice & 47 & , pssq_ice, pcd_du_ice, psen_ice & 48 & , pevp_ice, pwndm_ice, pfrac_oce & 49 & , ptaui_ice, ptauj_ice & 50 #endif 51 & ) 51 52 !--------------------------------------------------------------------------------------------------- 52 53 … … 54 55 !! *** ROUTINE abl_stp *** 55 56 !! 56 !! ** Purpose : Time-integration of the ABL model 57 !! ** Purpose : Time-integration of the ABL model 57 58 !! 58 !! ** Method : Compute atmospheric variables : vertical turbulence 59 !! ** Method : Compute atmospheric variables : vertical turbulence 59 60 !! + Coriolis term + newtonian relaxation 60 !! 61 !! 61 62 !! ** Action : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh 62 63 !! - Advance tracers to time n+1 (Euler backward scheme) … … 70 71 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psst ! sea-surface temperature [Celsius] 71 72 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu ! sea-surface u (U-point) 72 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv ! sea-surface v (V-point) 73 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv ! sea-surface v (V-point) 73 74 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq ! sea-surface humidity 74 75 REAL(wp) , INTENT(in ), DIMENSION(:,:,:) :: pu_dta ! large-scale windi … … 82 83 REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: psen ! Ch x Du 83 84 REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pevp ! Ce x Du 84 REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! ||uwnd|| 85 REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pwndm ! ||uwnd|| 85 86 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui ! taux 86 87 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj ! tauy 87 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaum ! ||tau|| 88 ! 89 #if defined key_si3 88 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaum ! ||tau|| 89 ! 90 #if defined key_si3 90 91 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptm_su ! ice-surface temperature [K] 91 92 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssu_ice ! ice-surface u (U-point) 92 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv_ice ! ice-surface v (V-point) 93 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq_ice ! ice-surface humidity 93 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssv_ice ! ice-surface v (V-point) 94 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pssq_ice ! ice-surface humidity 94 95 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pcd_du_ice ! Cd x Du over ice (T-point) 95 96 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: psen_ice ! Ch x Du over ice (T-point) 96 97 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pevp_ice ! Ce x Du over ice (T-point) 97 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndm_ice ! ||uwnd - uice|| 98 !REAL(wp) , INTENT(inout), DIMENSION(:,: ) :: pfrac_oce !!GS: out useless ? 99 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pfrac_oce ! 98 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndm_ice ! ||uwnd - uice|| 99 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pfrac_oce ! ocean fraction 100 100 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptaui_ice ! ice-surface taux stress (U-point) 101 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point) 102 #endif 103 ! 104 REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j 105 REAL(wp), DIMENSION(1:jpi,2:jpka ) :: zCF 106 REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: z_cft !--FL--to be removed after the test phase 107 ! 108 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a 109 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b 110 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_c 101 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point) 102 #endif 103 ! 104 REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j 105 REAL(wp), DIMENSION(1:jpi ,2:jpka) :: zCF 106 ! 107 REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_a 108 REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_b 109 REAL(wp), DIMENSION(1:jpi ,1:jpka) :: z_elem_c 111 110 ! 112 111 INTEGER :: ji, jj, jk, jtra, jbak ! dummy loop indices 113 112 REAL(wp) :: zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2 114 113 REAL(wp) :: zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce 115 ! 116 !!--------------------------------------------------------------------- 114 LOGICAL :: SemiImp_Cor = .TRUE. 115 ! 116 !!--------------------------------------------------------------------- 117 117 ! 118 118 IF(lwp .AND. kt == nit000) THEN ! control print … … 120 120 WRITE(numout,*) 'abl_stp : ABL time stepping' 121 121 WRITE(numout,*) '~~~~~~' 122 ENDIF 122 ENDIF 123 123 ! 124 124 IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) ) 125 !! Compute ustar squared as Cd || Uatm-Uoce ||^2 126 !! needed for surface boundary condition of TKE 125 IF( kt == nit000 ) ALLOCATE ( zrough( 1:jpi, 1:jpj ) ) 126 !! Compute ustar squared as Cd || Uatm-Uoce ||^2 127 !! needed for surface boundary condition of TKE 127 128 !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk) 128 129 DO_2D_11_11 129 130 zzoce = pCd_du (ji,jj) * pwndm (ji,jj) 130 131 #if defined key_si3 131 zzice = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 132 ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 132 zzice = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 133 ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 133 134 #else 134 ustar2(ji,jj) = zzoce 135 ustar2(ji,jj) = zzoce 135 136 #endif 137 zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(ji,jj), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce 136 138 END_2D 137 139 ! … … 140 142 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 141 143 142 CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj) 143 144 CALL abl_zdf_tke( ) !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj) 145 144 146 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 145 147 ! ! 2 *** Advance tracers to time n+1 146 148 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 147 149 148 150 !------------- 149 151 DO jj = 1, jpj ! outer loop !--> tq_abl computed on (1:jpi) x (1:jpj) 150 !------------- 151 ! Compute matrix elements for interior points 152 !------------- 153 ! Compute matrix elements for interior points 152 154 DO jk = 3, jpkam1 153 155 DO ji = 1, jpi ! vector opt. 154 z_elem_a( ji, jk) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal155 z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal156 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal157 END DO 158 END DO 159 ! Boundary conditions 160 DO ji = 1, jpi ! vector opt. 161 ! Neumann at the bottom 162 z_elem_a( ji, 2) = 0._wp163 z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 )156 z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 157 z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 158 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 159 END DO 160 END DO 161 ! Boundary conditions 162 DO ji = 1, jpi ! vector opt. 163 ! Neumann at the bottom 164 z_elem_a( ji, 2 ) = 0._wp 165 z_elem_c( ji, 2 ) = - rDt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 ) 164 166 ! Homogeneous Neumann at the top 165 z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )166 z_elem_c( ji, jpka) = 0._wp167 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji,jpka )168 END DO 167 z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 168 z_elem_c( ji, jpka ) = 0._wp 169 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 170 END DO 169 171 170 172 DO jtra = 1,jptq ! loop on active tracers 171 173 172 174 DO jk = 3, jpkam1 173 DO ji = 1,jpi 174 tq_abl ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl ( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side 175 !DO ji = 2, jpim1 176 DO ji = 1,jpi !!GS: to be checked if needed 177 tq_abl( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl( ji, jj, jk, nt_n, jtra ) ! initialize right-hand-side 175 178 END DO 176 179 END DO 177 180 178 181 IF(jtra == jp_ta) THEN 179 DO ji = 1,jpi ! boundary conditions for temperature180 zztmp1 = psen(ji, jj) 181 zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 ) 182 #if defined key_si3 183 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 182 DO ji = 1,jpi ! surface boundary condition for temperature 183 zztmp1 = psen(ji, jj) 184 zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 ) 185 #if defined key_si3 186 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 184 187 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj) 185 #endif 186 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1187 tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2 , nt_n, jtra ) + rDt_abl * zztmp2188 #endif 189 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 190 tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2 188 191 tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra ) 189 END DO 190 ELSE 191 DO ji = 1,jpi ! boundary conditions for humidity192 zztmp1 = pevp(ji, jj) 193 zztmp2 = pevp(ji, jj) * pssq(ji, jj) 194 #if defined key_si3 195 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 196 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj) 197 #endif 198 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2) + rDt_abl * zztmp1199 tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2 , nt_n, jtra ) + rDt_abl * zztmp2192 END DO 193 ELSE ! jp_qa 194 DO ji = 1,jpi ! surface boundary condition for humidity 195 zztmp1 = pevp(ji, jj) 196 zztmp2 = pevp(ji, jj) * pssq(ji, jj) 197 #if defined key_si3 198 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 199 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj) 200 #endif 201 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 202 tq_abl ( ji, jj, 2 , nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rDt_abl * zztmp2 200 203 tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra ) 201 END DO 204 END DO 202 205 END IF 203 206 !! 204 207 !! Matrix inversion 205 208 !! ---------------------------------------------------------- 206 DO ji = 1,jpi 207 zcff = 1._wp / z_elem_b( ji, 2 )208 zCF ( ji, 2) = - zcff * z_elem_c( ji, 2 )209 tq_abl( ji,jj,2,nt_a,jtra) = zcff * tq_abl(ji,jj,2,nt_a,jtra)210 END DO 211 212 DO jk = 3, jpka 213 DO ji = 1,jpi 214 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) )209 DO ji = 1,jpi 210 zcff = 1._wp / z_elem_b( ji, 2 ) 211 zCF ( ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 212 tq_abl( ji, jj, 2, nt_a, jtra ) = zcff * tq_abl( ji, jj, 2, nt_a, jtra ) 213 END DO 214 215 DO jk = 3, jpka 216 DO ji = 1,jpi 217 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF( ji, jk-1 ) ) 215 218 zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) 216 219 tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk ,nt_a,jtra) & 217 & - z_elem_a(ji, jk) *tq_abl(ji,jj,jk-1,nt_a,jtra) )218 END DO 219 END DO 220 !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ... 221 DO jk = jpkam1,2,-1 222 DO ji = 1,jpi 220 & - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) ) 221 END DO 222 END DO 223 !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ... 224 DO jk = jpkam1,2,-1 225 DO ji = 1,jpi 223 226 tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) + & 224 227 & zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra) 225 228 END DO 226 229 END DO 227 228 END DO !<-- loop on tracers 229 !! 230 !------------- 231 END DO ! end outer loop 232 !------------- 233 234 230 231 END DO !<-- loop on tracers 232 !! 233 !------------- 234 END DO ! end outer loop 235 !------------- 236 235 237 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 236 238 ! ! 3 *** Compute Coriolis term with geostrophic guide 237 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 238 !------------- 239 DO jk = 2, jpka ! outer loop 240 !------------- 241 ! 242 ! Advance u_abl & v_abl to time n+1 243 DO_2D_11_11 244 zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl ) ! (f dt)**2 239 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 240 IF( SemiImp_Cor ) THEN 241 242 !------------- 243 DO jk = 2, jpka ! outer loop 244 !------------- 245 ! 246 ! Advance u_abl & v_abl to time n+1 247 DO_2D_11_11 248 zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl ) ! (f dt)**2 249 250 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 251 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * u_abl( ji, jj, jk, nt_n ) & 252 & + rDt_abl * fft_abl(ji, jj) * v_abl( ji, jj, jk, nt_n ) ) & 253 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 245 254 246 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 247 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n ) & 248 & + rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj , jk, nt_n ) ) & 249 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 250 251 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 252 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n ) & 253 & - rDt_abl * fft_abl(ji, jj) * u_abl ( ji , jj, jk, nt_n ) ) & 254 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 255 END_2D 256 ! 257 !------------- 258 END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 259 !------------- 260 ! 261 IF( ln_geos_winds ) THEN 262 DO jj = 1, jpj ! outer loop 263 DO jk = 1, jpka 264 DO ji = 1, jpi 265 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) & 266 & - rDt_abl * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk) 267 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) & 268 & + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj ) * pgu_dta(ji ,jj ,jk) 269 END DO 270 END DO 271 END DO 272 END IF 273 !------------- 274 ! 275 IF( ln_hpgls_frc ) THEN 276 DO jj = 1, jpj ! outer loop 277 DO jk = 1, jpka 278 DO ji = 1, jpi 279 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 280 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 255 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 256 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * v_abl( ji, jj, jk, nt_n ) & 257 & - rDt_abl * fft_abl(ji, jj) * u_abl( ji, jj, jk, nt_n ) ) & 258 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 259 END_2D 260 ! 261 !------------- 262 END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 263 !------------- 264 ! 265 IF( ln_geos_winds ) THEN 266 DO jj = 1, jpj ! outer loop 267 DO jk = 1, jpka 268 DO ji = 1, jpi 269 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) & 270 & - rDt_abl * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk) 271 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) & 272 & + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj ) * pgu_dta(ji ,jj ,jk) 273 END DO 274 END DO 275 END DO 276 END IF 277 ! 278 IF( ln_hpgls_frc ) THEN 279 DO jj = 1, jpj ! outer loop 280 DO jk = 1, jpka 281 DO ji = 1, jpi 282 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 283 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 284 ENDDO 281 285 ENDDO 282 286 ENDDO 283 ENDDO 284 END IF 285 286 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 287 ! ! 4 *** Advance u,v to time n+1 288 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 289 ! 290 ! Vertical diffusion for u_abl 287 END IF 288 289 ELSE ! SemiImp_Cor = .FALSE. 290 291 IF( ln_geos_winds ) THEN 292 293 !------------- 294 DO jk = 2, jpka ! outer loop 295 !------------- 296 ! 297 IF( MOD( kt, 2 ) == 0 ) then 298 ! Advance u_abl & v_abl to time n+1 299 DO jj = 1, jpj 300 DO ji = 1, jpi 301 zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_n ) - pgv_dta(ji ,jj ,jk) ) 302 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff 303 zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_a ) - pgu_dta(ji ,jj ,jk) ) 304 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff ) 305 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * u_abl( ji, jj, jk, nt_a ) 306 END DO 307 END DO 308 ELSE 309 DO jj = 1, jpj 310 DO ji = 1, jpi 311 zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_n ) - pgu_dta(ji ,jj ,jk) ) 312 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff 313 zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_a ) - pgv_dta(ji ,jj ,jk) ) 314 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff ) 315 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * v_abl( ji, jj, jk, nt_a ) 316 END DO 317 END DO 318 END IF 319 ! 320 !------------- 321 END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 322 !------------- 323 324 ENDIF ! ln_geos_winds 325 326 ENDIF ! SemiImp_Cor 327 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 328 ! ! 4 *** Advance u,v to time n+1 329 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 330 ! 331 ! Vertical diffusion for u_abl 291 332 !------------- 292 333 DO jj = 1, jpj ! outer loop 293 !------------- 334 !------------- 294 335 295 336 DO jk = 3, jpkam1 296 DO ji = 1, jpi 297 z_elem_a( ji, 298 z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal299 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )! diagonal300 END DO 301 END DO 302 303 DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi) 337 DO ji = 1, jpi 338 z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 339 z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 340 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 341 END DO 342 END DO 343 344 DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi) 304 345 !++ Surface boundary condition 305 z_elem_a( ji, 2) = 0._wp306 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )307 ! 308 309 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) 310 #if defined key_si3 346 z_elem_a( ji, 2 ) = 0._wp 347 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 348 ! 349 zztmp1 = pcd_du(ji, jj) 350 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) * rn_vfac 351 #if defined key_si3 311 352 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) 312 zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) ) 313 314 #endif 315 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 353 zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj) ) * rn_vfac 354 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 355 #endif 356 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 316 357 u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2 317 318 !++ Top Neumann B.C. 319 !z_elem_a( ji, jpka ) = - 0.5_wp * rDt_abl * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka ) 320 !z_elem_c( ji, jpka ) = 0._wp 321 !z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 322 !++ Top Dirichlet B.C. 323 z_elem_a( ji, jpka ) = 0._wp 324 z_elem_c( ji, jpka ) = 0._wp 325 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 326 u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 327 END DO 358 359 ! idealized test cases only 360 !IF( ln_topbc_neumann ) THEN 361 ! !++ Top Neumann B.C. 362 ! z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 363 ! z_elem_c( ji, jpka ) = 0._wp 364 ! z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 365 ! !u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl ( ji, jj, jpka, nt_a ) 366 !ELSE 367 !++ Top Dirichlet B.C. 368 z_elem_a( ji, jpka ) = 0._wp 369 z_elem_c( ji, jpka ) = 0._wp 370 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 371 u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 372 !ENDIF 373 374 END DO 328 375 !! 329 376 !! Matrix inversion 330 377 !! ---------------------------------------------------------- 331 DO ji = 2, jpi 378 !DO ji = 2, jpi 379 DO ji = 1, jpi !!GS: TBI 332 380 zcff = 1._wp / z_elem_b( ji, 2 ) 333 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 381 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 334 382 u_abl (ji,jj,2,nt_a) = zcff * u_abl(ji,jj,2,nt_a) 335 END DO 336 337 DO jk = 3, jpka 338 DO ji = 2, jpi 383 END DO 384 385 DO jk = 3, jpka 386 DO ji = 2, jpi 339 387 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) ) 340 388 zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) … … 343 391 END DO 344 392 END DO 345 346 DO jk = jpkam1,2,-1 393 394 DO jk = jpkam1,2,-1 347 395 DO ji = 2, jpi 348 396 u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a) 349 397 END DO 350 398 END DO 351 352 !------------- 353 END DO ! end outer loop 354 !------------- 355 356 ! 357 ! Vertical diffusion for v_abl 399 400 !------------- 401 END DO ! end outer loop 402 !------------- 403 404 ! 405 ! Vertical diffusion for v_abl 358 406 !------------- 359 407 DO jj = 2, jpj ! outer loop 360 !------------- 408 !------------- 361 409 ! 362 410 DO jk = 3, jpkam1 363 DO ji = 1, jpi 364 z_elem_a( ji, 365 z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal366 z_elem_b( ji, 367 END DO 368 END DO 369 370 DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 411 DO ji = 1, jpi 412 z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 413 z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 414 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 415 END DO 416 END DO 417 418 DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 371 419 !++ Surface boundary condition 372 z_elem_a( ji, 2) = 0._wp373 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )374 ! 375 376 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) 377 #if defined key_si3 420 z_elem_a( ji, 2 ) = 0._wp 421 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 422 ! 423 zztmp1 = pcd_du(ji, jj) 424 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) * rn_vfac 425 #if defined key_si3 378 426 zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) 379 zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) ) 380 381 #endif 382 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 427 zzice = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) ) * rn_vfac 428 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 429 #endif 430 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 383 431 v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2 384 !++ Top Neumann B.C. 385 !z_elem_a( ji, jpka ) = -rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 386 !z_elem_c( ji, jpka ) = 0._wp 387 !z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 388 !++ Top Dirichlet B.C. 389 z_elem_a( ji, jpka ) = 0._wp 390 z_elem_c( ji, jpka ) = 0._wp 391 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 392 v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 393 END DO 432 433 ! idealized test cases only 434 !IF( ln_topbc_neumann ) THEN 435 ! !++ Top Neumann B.C. 436 ! z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 437 ! z_elem_c( ji, jpka ) = 0._wp 438 ! z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 439 ! !v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl ( ji, jj, jpka, nt_a ) 440 !ELSE 441 !++ Top Dirichlet B.C. 442 z_elem_a( ji, jpka ) = 0._wp 443 z_elem_c( ji, jpka ) = 0._wp 444 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 445 v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 446 !ENDIF 447 448 END DO 394 449 !! 395 450 !! Matrix inversion 396 451 !! ---------------------------------------------------------- 397 DO ji = 1, jpi 452 DO ji = 1, jpi 398 453 zcff = 1._wp / z_elem_b( ji, 2 ) 399 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 400 v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a ) 401 END DO 402 403 DO jk = 3, jpka 404 DO ji = 1, jpi 454 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 455 v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a ) 456 END DO 457 458 DO jk = 3, jpka 459 DO ji = 1, jpi 405 460 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) ) 406 461 zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) … … 409 464 END DO 410 465 END DO 411 412 DO jk = jpkam1,2,-1 413 DO ji = 1, jpi 466 467 DO jk = jpkam1,2,-1 468 DO ji = 1, jpi 414 469 v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a) 415 470 END DO 416 471 END DO 417 ! 418 !------------- 419 END DO ! end outer loop 420 !------------- 421 422 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 423 ! ! 5 *** Apply nudging on the dynamics and the tracers 424 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 425 z_cft(:,:,:) = 0._wp 426 472 ! 473 !------------- 474 END DO ! end outer loop 475 !------------- 476 477 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 478 ! ! 5 *** Apply nudging on the dynamics and the tracers 479 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 480 427 481 IF( nn_dyn_restore > 0 ) THEN 428 !------------- 482 !------------- 429 483 DO jk = 2, jpka ! outer loop 430 !------------- 484 !------------- 431 485 DO_2D_01_01 432 486 zcff1 = pblh( ji, jj ) 433 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 434 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 487 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 488 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 435 489 zmsk = msk_abl(ji,jj) 436 490 zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2 & 437 491 & + jp_alp1_dyn * zsig + jp_alp0_dyn 438 492 zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points 439 ! rn_Dt = rDt_abl / nn_fsbc 493 ! rn_Dt = rDt_abl / nn_fsbc 440 494 zcff = zcff * rest_eq(ji,jj) 441 z_cft( ji, jj, jk ) = zcff442 495 u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) & 443 & + zcff * pu_dta( ji, jj, jk ) 496 & + zcff * pu_dta( ji, jj, jk ) 444 497 v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * v_abl( ji, jj, jk, nt_a ) & 445 498 & + zcff * pv_dta( ji, jj, jk ) … … 447 500 !------------- 448 501 END DO ! end outer loop 449 !------------- 502 !------------- 450 503 END IF 451 504 452 !------------- 505 !------------- 453 506 DO jk = 2, jpka ! outer loop 454 !------------- 507 !------------- 455 508 DO_2D_11_11 456 509 zcff1 = pblh( ji, jj ) 457 510 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 458 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 511 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 459 512 zmsk = msk_abl(ji,jj) 460 513 zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2 & 461 514 & + jp_alp1_tra * zsig + jp_alp0_tra 462 515 zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points 463 ! rn_Dt = rDt_abl / nn_fsbc 464 !z_cft( ji, jj, jk ) = zcff 516 ! rn_Dt = rDt_abl / nn_fsbc 465 517 tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) & 466 518 & + zcff * pt_dta( ji, jj, jk ) 467 519 468 520 tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa ) & 469 521 & + zcff * pq_dta( ji, jj, jk ) 470 522 471 523 END_2D 472 524 !------------- 473 525 END DO ! end outer loop 474 !------------- 475 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 476 ! ! 6 *** MPI exchanges 477 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 478 ! 479 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1., v_abl(:,:,:,nt_a ), 'T', -1.)526 !------------- 527 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 528 ! ! 6 *** MPI exchanges 529 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 530 ! 531 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1., v_abl(:,:,:,nt_a) , 'T', -1. ) 480 532 CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T', 1., kfillmode = jpfillnothing ) ! ++++ this should not be needed... 481 533 ! 482 ! first ABL level 534 #if defined key_iomput 535 ! 2D & first ABL level 536 IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh", pblh(:,: ) ) 483 537 IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl", u_abl(:,:,2,nt_a ) ) 484 538 IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl", v_abl(:,:,2,nt_a ) ) … … 489 543 IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta", pt_dta(:,:,2 ) ) 490 544 IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta", pq_dta(:,:,2 ) ) 491 ! all ABL levels 492 IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) ) 493 IF ( iom_use("v_abl" ) ) CALL iom_put ( "v_abl" , v_abl(:,:,2:jpka,nt_a ) ) 494 IF ( iom_use("t_abl" ) ) CALL iom_put ( "t_abl" , tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 495 IF ( iom_use("q_abl" ) ) CALL iom_put ( "q_abl" , tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 496 IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a ) ) 497 IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka ) ) 498 IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka ) ) 499 IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxl_abl(:,:,2:jpka ) ) 500 IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh" , pblh(:,: ) ) 501 ! debug (to be removed) 545 ! debug 2D 546 IF( ln_geos_winds ) THEN 547 IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) ) 548 IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) ) 549 END IF 550 IF( ln_hpgls_frc ) THEN 551 IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 552 IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 553 END IF 554 ! 3D (all ABL levels) 555 IF ( iom_use("u_abl" ) ) CALL iom_put ( "u_abl" , u_abl(:,:,2:jpka,nt_a ) ) 556 IF ( iom_use("v_abl" ) ) CALL iom_put ( "v_abl" , v_abl(:,:,2:jpka,nt_a ) ) 557 IF ( iom_use("t_abl" ) ) CALL iom_put ( "t_abl" , tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 558 IF ( iom_use("q_abl" ) ) CALL iom_put ( "q_abl" , tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 559 IF ( iom_use("tke_abl" ) ) CALL iom_put ( "tke_abl" , tke_abl(:,:,2:jpka,nt_a ) ) 560 IF ( iom_use("avm_abl" ) ) CALL iom_put ( "avm_abl" , avm_abl(:,:,2:jpka ) ) 561 IF ( iom_use("avt_abl" ) ) CALL iom_put ( "avt_abl" , avt_abl(:,:,2:jpka ) ) 562 IF ( iom_use("mxlm_abl") ) CALL iom_put ( "mxlm_abl", mxlm_abl(:,:,2:jpka ) ) 563 IF ( iom_use("mxld_abl") ) CALL iom_put ( "mxld_abl", mxld_abl(:,:,2:jpka ) ) 564 ! debug 3D 502 565 IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta", pu_dta(:,:,2:jpka) ) 503 566 IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta", pv_dta(:,:,2:jpka) ) 504 567 IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta", pt_dta(:,:,2:jpka) ) 505 568 IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta", pq_dta(:,:,2:jpka) ) 506 IF ( iom_use("coeft") ) CALL iom_put ( "coeft", z_cft(:,:,2:jpka) )507 569 IF( ln_geos_winds ) THEN 508 IF ( iom_use("u z1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) )509 IF ( iom_use("v z1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) )570 IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka) ) 571 IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgv_dta(:,:,2:jpka) ) 510 572 END IF 511 573 IF( ln_hpgls_frc ) THEN 512 IF ( iom_use("u z1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp))513 IF ( iom_use("v z1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp))574 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) ) 575 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) ) 514 576 END IF 515 ! 577 #endif 516 578 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 517 579 ! ! 7 *** Finalize flux computation 518 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 519 580 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 581 ! 520 582 DO_2D_11_11 521 ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta ) 522 zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa ) 523 !zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ... 524 ! & ( R_dry*ztemp * ( 1._wp + rctv0*zhumi ) ) 525 zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 526 psen ( ji, jj ) = cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) 527 pevp ( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) ) 528 rhoa( ji, jj ) = zcff 583 ztemp = tq_abl( ji, jj, 2, nt_a, jp_ta ) 584 zhumi = tq_abl( ji, jj, 2, nt_a, jp_qa ) 585 zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 586 psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) !GS: negative sign to respect aerobulk convention 587 pevp( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) ) 588 rhoa( ji, jj ) = zcff 529 589 END_2D 530 590 531 591 DO_2D_01_01 592 <<<<<<< .working 532 593 zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * ( pssu(ji ,jj) + pssu(ji-1,jj) ) 533 594 zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * ( pssv(ji,jj ) + pssv(ji,jj-1) ) 595 ||||||| .merge-left.r13136 596 zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji ,jj) + pssu(ji-1,jj) ) 597 zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj ) + pssv(ji,jj-1) ) 598 ======= 599 zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji,jj) + pssu(ji-1,jj ) ) 600 zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj) + pssv(ji ,jj-1) ) 601 >>>>>>> .merge-right.r13211 534 602 END_2D 535 ! 603 ! 536 604 CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. ) 537 605 ! … … 539 607 DO_2D_11_11 540 608 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 541 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) )! * msk_abl(ji,jj)609 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 542 610 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 543 611 544 612 pwndm (ji,jj) = zcff 545 613 ptaum (ji,jj) = zztmp * zcff … … 564 632 565 633 IF(sn_cfctl%l_prtctl) THEN 566 CALL prt_ctl( tab2d_1=p wndm , clinfo1=' abl_stp: wndm : ' )567 CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ')568 CALL prt_ctl( tab2d_ 2=ptauj , clinfo2= 'vtau: ' )634 CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ', mask1=umask, & 635 & tab2d_2=ptauj , clinfo2=' vtau : ', mask2=vmask ) 636 CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' ) 569 637 ENDIF 570 638 571 639 #if defined key_si3 640 <<<<<<< .working 572 641 ! ------------------------------------------------------------ ! 573 642 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! … … 583 652 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' & 584 653 & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' ) 654 ||||||| .merge-left.r13136 655 ! ------------------------------------------------------------ ! 656 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 657 ! ------------------------------------------------------------ ! 658 DO_2D_00_00 659 660 zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 661 zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 662 663 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 664 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 665 & * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 666 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 667 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 668 & * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 669 END_2D 670 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 671 ! 672 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' & 673 & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' ) 674 ======= 675 ! ------------------------------------------------------------ ! 676 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 677 ! ------------------------------------------------------------ ! 678 DO_2D_00_00 679 680 zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 681 zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 682 683 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 684 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 685 & * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 686 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 687 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 688 & * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 689 END_2D 690 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 691 ! 692 IF(sn_cfctl%l_prtctl) THEN 693 CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask, & 694 & tab2d_2=ptauj_ice , clinfo2=' vtau_ice : ', mask2=vmask ) 695 END IF 696 >>>>>>> .merge-right.r13211 585 697 #endif 586 698 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 615 727 !! (= Kz dz[Ub] * dz[Un] ) 616 728 !! --------------------------------------------------------------------- 617 INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn 729 INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn 618 730 INTEGER, DIMENSION(1:jpi ) :: ikbl 619 731 REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv 620 REAL(wp) :: zdU , zdV, zcff1,zshear,zbuoy,zsig, zustar2621 REAL(wp) :: zdU2, zdV2622 REAL(wp) :: zwndi, zwndj732 REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2 733 REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 ! zbuoy for BL89 734 REAL(wp) :: zwndi, zwndj 623 735 REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2 624 736 REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2 625 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF 737 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF 626 738 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a 627 739 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b … … 629 741 LOGICAL :: ln_Patankar = .FALSE. 630 742 LOGICAL :: ln_dumpvar = .FALSE. 631 LOGICAL , DIMENSION(1:jpi ) :: ln_foundl 743 LOGICAL , DIMENSION(1:jpi ) :: ln_foundl 632 744 ! 633 745 tind = nt_n … … 641 753 !------------- 642 754 ! 643 ! Compute vertical shear 755 ! Compute vertical shear 644 756 DO jk = 2, jpkam1 645 DO ji = 1, jpi646 zcff = 1.0_wp / e3w_abl( jk )**2 647 zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 757 DO ji = 1, jpi 758 zcff = 1.0_wp / e3w_abl( jk )**2 759 zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 648 760 zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 649 zsh2(ji,jk) = zdU+zdV 650 END DO 651 END DO 761 zsh2(ji,jk) = zdU+zdV !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 ) 762 END DO 763 END DO 652 764 ! 653 765 ! Compute brunt-vaisala frequency 654 766 DO jk = 2, jpkam1 655 DO ji = 1,jpi 656 zcff = grav * itvref / e3w_abl( jk ) 767 DO ji = 1,jpi 768 zcff = grav * itvref / e3w_abl( jk ) 657 769 zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta) 658 770 zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) & … … 660 772 zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi) 661 773 END DO 662 END DO 774 END DO 663 775 ! 664 776 ! Terms for the tridiagonal problem 665 777 DO jk = 2, jpkam1 666 DO ji = 1, jpi667 zshear = zsh2( ji, jk )! zsh2 is already multiplied by Avm_abl at this point668 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation669 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )670 671 z_elem_a( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal672 z_elem_c( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal778 DO ji = 1, jpi 779 zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point 780 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 781 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 782 783 z_elem_a( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk ) + Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal 784 z_elem_c( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk ) + Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal 673 785 IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE 674 z_elem_b( ji, jk )= e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &675 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) ! diagonal676 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )! right-hand-side786 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 787 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) ! diagonal 788 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) ) ! right-hand-side 677 789 ELSE 678 z_elem_b( ji, jk )= e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) &679 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) & ! diagonal680 & - e3w_abl(jk) * rDt_abl * zbuoy681 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side790 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 791 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) & ! diagonal 792 & - e3w_abl(jk) * rDt_abl * zbuoy 793 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side 682 794 END IF 683 795 END DO 684 END DO 685 686 DO ji = 1,jpi ! vector opt. 687 zesrf = MAX( 4.63_wp * ustar2(ji,jj), tke_min ) 688 zetop = tke_min 689 z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp 690 z_elem_a ( ji, jpka ) = 0._wp; z_elem_c ( ji, jpka ) = 0._wp; z_elem_b ( ji, jpka ) = 1._wp 691 tke_abl( ji, jj, 1, nt_a ) = zesrf 692 tke_abl( ji, jj, jpka, nt_a ) = zetop 693 zbn2(ji,jj, 1) = zbn2( ji,jj, 2) 694 zsh2(ji, 1) = zsh2( ji, 2) 695 zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 696 zsh2(ji, jpka) = zsh2( ji , jpkam1) 697 END DO 796 END DO 797 798 DO ji = 1,jpi ! vector opt. 799 zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min ) 800 zetop = tke_min 801 802 z_elem_a ( ji, 1 ) = 0._wp 803 z_elem_c ( ji, 1 ) = 0._wp 804 z_elem_b ( ji, 1 ) = 1._wp 805 tke_abl ( ji, jj, 1, nt_a ) = zesrf 806 807 !++ Top Neumann B.C. 808 !z_elem_a ( ji, jpka ) = - 0.5 * rDt_abl * rn_Sch * (Avm_abl(ji,jj, jpka-1 )+Avm_abl(ji,jj, jpka )) / e3t_abl( jpka ) 809 !z_elem_c ( ji, jpka ) = 0._wp 810 !z_elem_b ( ji, jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 811 !tke_abl ( ji, jj, jpka, nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 812 813 !++ Top Dirichlet B.C. 814 z_elem_a ( ji, jpka ) = 0._wp 815 z_elem_c ( ji, jpka ) = 0._wp 816 z_elem_b ( ji, jpka ) = 1._wp 817 tke_abl ( ji, jj, jpka, nt_a ) = zetop 818 819 zbn2 ( ji, jj, 1 ) = zbn2 ( ji, jj, 2 ) 820 zsh2 ( ji, 1 ) = zsh2 ( ji, 2 ) 821 zbn2 ( ji, jj, jpka ) = zbn2 ( ji, jj, jpkam1 ) 822 zsh2 ( ji, jpka ) = zsh2 ( ji , jpkam1 ) 823 END DO 698 824 !! 699 825 !! Matrix inversion … … 701 827 DO ji = 1,jpi 702 828 zcff = 1._wp / z_elem_b( ji, 1 ) 703 zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 ) 704 tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a ) 705 END DO 706 707 DO jk = 2, jpka 829 zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 ) 830 tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a ) 831 END DO 832 833 DO jk = 2, jpka 708 834 DO ji = 1,jpi 709 835 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) ) … … 713 839 END DO 714 840 END DO 715 716 DO jk = jpkam1,1,-1 841 842 DO jk = jpkam1,1,-1 717 843 DO ji = 1,jpi 718 844 tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a) 719 845 END DO 720 846 END DO 721 722 !!FL should not be needed because of Patankar procedure 847 848 !!FL should not be needed because of Patankar procedure 723 849 tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min ) 724 850 … … 726 852 !! Diagnose PBL height 727 853 !! ---------------------------------------------------------- 728 729 730 ! 854 855 856 ! 731 857 ! arrays zRH, zFC and zCF are available at this point 732 858 ! and zFC(:, 1 ) = 0. 733 859 ! diagnose PBL height based on zsh2 and zbn2 734 860 zFC ( : ,1) = 0._wp 735 ikbl( 1:jpi ) = 0 736 861 ikbl( 1:jpi ) = 0 862 737 863 DO jk = 2,jpka 738 DO ji = 1, jpi 864 DO ji = 1, jpi 739 865 zcff = ghw_abl( jk-1 ) 740 866 zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) ) … … 762 888 ELSE 763 889 pblh( ji, jj ) = ghw_abl(jpka) 764 END IF 765 END DO 766 !------------- 767 END DO 768 !------------- 769 ! 770 ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 890 END IF 891 END DO 892 !------------- 893 END DO 894 !------------- 895 ! 896 ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 771 897 IF(ln_smth_pblh) THEN 772 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) 898 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing) 773 899 CALL smooth_pblh( pblh, msk_abl ) 774 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) 900 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing) 775 901 ENDIF 776 902 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 780 906 SELECT CASE ( nn_amxl ) 781 907 ! 782 CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom 783 # define zlup zRH 784 # define zldw zFC 908 CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom 909 # define zlup zRH 910 # define zldw zFC 785 911 DO jj = 1, jpj ! outer loop 786 912 ! 787 913 DO ji = 1, jpi 788 mxl_abl ( ji, jj, 1 ) = 0._wp 789 mxl_abl ( ji, jj, jpka ) = mxl_min 790 zldw( ji, 1 ) = 0._wp 791 zlup( ji, jpka ) = 0._wp 792 END DO 793 ! 794 DO jk = 2, jpkam1 795 DO ji = 1, jpi 796 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 797 mxl_abl( ji, jj, jk ) = MAX( mxl_min, & 798 & SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) ) 799 END DO 800 END DO 914 mxld_abl( ji, jj, 1 ) = mxl_min 915 mxld_abl( ji, jj, jpka ) = mxl_min 916 mxlm_abl( ji, jj, 1 ) = mxl_min 917 mxlm_abl( ji, jj, jpka ) = mxl_min 918 zldw ( ji, 1 ) = zrough(ji,jj) * rn_Lsfc 919 zlup ( ji, jpka ) = mxl_min 920 END DO 921 ! 922 DO jk = 2, jpkam1 923 DO ji = 1, jpi 924 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 925 mxlm_abl( ji, jj, jk ) = MAX( mxl_min, & 926 & SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) ) 927 END DO 928 END DO 801 929 ! 802 930 ! Limit mxl 803 DO jk = jpkam1,1,-1 804 DO ji = 1, jpi 805 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl _abl(ji, jj, jk) )806 END DO 807 END DO 931 DO jk = jpkam1,1,-1 932 DO ji = 1, jpi 933 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 934 END DO 935 END DO 808 936 ! 809 937 DO jk = 2, jpka 810 DO ji = 1, jpi 811 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) ) 812 END DO 813 END DO 938 DO ji = 1, jpi 939 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 940 END DO 941 END DO 942 ! 943 ! DO jk = 1, jpka 944 ! DO ji = 1, jpi 945 ! mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 946 ! mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ), zlup( ji, jk ) ) 947 ! END DO 948 ! END DO 814 949 ! 815 950 DO jk = 1, jpka 816 951 DO ji = 1, jpi 817 mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 818 END DO 819 END DO 820 ! 821 END DO 822 # undef zlup 823 # undef zldw 824 ! 825 ! 826 CASE ( 1 ) ! length-scale computed as the distance to the PBL height 827 DO jj = 1,jpj ! outer loop 828 ! 829 DO ji = 1, jpi ! vector opt. 830 zcff = 1._wp / pblh( ji, jj ) ! inverse of hbl 831 DO jk = 1, jpka 832 zsig = MIN( zcff * ghw_abl( jk ), 1. ) 833 zcff1 = pblh( ji, jj ) 834 mxl_abl( ji, jj, jk ) = mxl_min & 835 & + zsig * ( amx1*zcff1 + bmx1*mxl_min ) & 836 & + zsig * zsig * ( amx2*zcff1 + bmx2*mxl_min ) & 837 & + zsig**3 * ( amx3*zcff1 + bmx3*mxl_min ) & 838 & + zsig**4 * ( amx4*zcff1 + bmx4*mxl_min ) & 839 & + zsig**5 * ( amx5*zcff1 + bmx5*mxl_min ) 840 END DO 841 END DO 842 ! 843 END DO 952 ! zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 953 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 954 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 955 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 956 END DO 957 END DO 958 ! 959 END DO 960 # undef zlup 961 # undef zldw 962 ! 963 ! 964 CASE ( 1 ) ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom 965 # define zlup zRH 966 # define zldw zFC 967 DO jj = 1, jpj ! outer loop 968 ! 969 DO jk = 2, jpkam1 970 DO ji = 1,jpi 971 zcff = 1.0_wp / e3w_abl( jk )**2 972 zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 973 zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 974 zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 975 ENDDO 976 ENDDO 977 ! 978 DO ji = 1, jpi 979 zcff = zrough(ji,jj) * rn_Lsfc 980 mxld_abl ( ji, jj, 1 ) = zcff 981 mxld_abl ( ji, jj, jpka ) = mxl_min 982 mxlm_abl ( ji, jj, 1 ) = zcff 983 mxlm_abl ( ji, jj, jpka ) = mxl_min 984 zldw ( ji, 1 ) = zcff 985 zlup ( ji, jpka ) = mxl_min 986 END DO 987 ! 988 DO jk = 2, jpkam1 989 DO ji = 1, jpi 990 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 991 zcff = 2.*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) & 992 & + SQRT( rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.*zbuoy ) ) 993 mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff ) 994 END DO 995 END DO 996 ! 997 ! Limit mxl 998 DO jk = jpkam1,1,-1 999 DO ji = 1, jpi 1000 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 1001 END DO 1002 END DO 1003 ! 1004 DO jk = 2, jpka 1005 DO ji = 1, jpi 1006 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 1007 END DO 1008 END DO 1009 ! 1010 DO jk = 1, jpka 1011 DO ji = 1, jpi 1012 !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1013 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1014 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1015 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1016 !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) ) 1017 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1018 END DO 1019 END DO 1020 ! 1021 END DO 1022 # undef zlup 1023 # undef zldw 844 1024 ! 845 1025 CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale 846 1026 ! 847 # define zlup zRH 848 # define zldw zFC 1027 # define zlup zRH 1028 # define zldw zFC 849 1029 ! zCF is used for matrix inversion 850 ! 1030 ! 851 1031 DO jj = 1, jpj ! outer loop 852 853 DO ji = 1, jpi 854 zlup( ji, 1 ) = mxl_min 855 zldw( ji, 1 ) = mxl_min 1032 1033 DO ji = 1, jpi 1034 zcff = zrough(ji,jj) * rn_Lsfc 1035 zlup( ji, 1 ) = zcff 1036 zldw( ji, 1 ) = zcff 856 1037 zlup( ji, jpka ) = mxl_min 857 zldw( ji, jpka ) = mxl_min 858 END DO 859 1038 zldw( ji, jpka ) = mxl_min 1039 END DO 1040 860 1041 DO jk = 2,jpka-1 861 1042 DO ji = 1, jpi 862 1043 zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 863 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 864 END DO 865 END DO 1044 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 1045 END DO 1046 END DO 866 1047 !! 867 1048 !! BL89 search for lup 868 !! ---------------------------------------------------------- 869 DO jk=2,jpka-1 1049 !! ---------------------------------------------------------- 1050 DO jk=2,jpka-1 870 1051 ! 871 1052 DO ji = 1, jpi … … 873 1054 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 874 1055 ln_foundl(ji ) = .false. 875 END DO 876 ! 1056 END DO 1057 ! 877 1058 DO jkup=jk+1,jpka-1 878 1059 DO ji = 1, jpi 1060 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1061 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 879 1062 zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 880 & ( zb n2(ji,jj,jkup )*(ghw_abl(jkup )-ghw_abl(jk)) &881 & + zb n2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )1063 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & 1064 & + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) 882 1065 IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 883 1066 zcff2 = ghw_abl(jkup ) - ghw_abl(jk) 884 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 1067 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 885 1068 zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / & 886 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 887 zlup(ji,jk) = zcff 1069 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 1070 zlup(ji,jk) = zcff 1071 zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk) 888 1072 ln_foundl(ji) = .true. 889 1073 END IF … … 891 1075 END DO 892 1076 ! 893 END DO 1077 END DO 894 1078 !! 895 1079 !! BL89 search for ldwn 896 !! ---------------------------------------------------------- 897 DO jk=2,jpka-1 1080 !! ---------------------------------------------------------- 1081 DO jk=2,jpka-1 898 1082 ! 899 1083 DO ji = 1, jpi … … 901 1085 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 902 1086 ln_foundl(ji ) = .false. 903 END DO 904 ! 1087 END DO 1088 ! 905 1089 DO jkdwn=jk-1,1,-1 906 DO ji = 1, jpi 1090 DO ji = 1, jpi 1091 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1092 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 907 1093 zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) & 908 & * ( zb n2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &909 + zb n2(ji,jj,jkdwn )*(ghw_abl(jk)-ghw_abl(jkdwn )) )910 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1094 & * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 1095 + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) 1096 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 911 1097 zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 912 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 1098 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 913 1099 zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / & 914 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 915 zldw(ji,jk) = zcff 916 ln_foundl(ji) = .true. 917 END IF 918 END DO 919 END DO 920 ! 1100 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1101 zldw(ji,jk) = zcff 1102 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1103 ln_foundl(ji) = .true. 1104 END IF 1105 END DO 1106 END DO 1107 ! 921 1108 END DO 922 1109 923 1110 DO jk = 1, jpka 924 DO ji = 1, jpi 925 mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min ) 926 END DO 927 END DO 1111 DO ji = 1, jpi 1112 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1113 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1114 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1115 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1116 END DO 1117 END DO 928 1118 929 1119 END DO 930 # undef zlup 931 # undef zldw 932 ! 933 END SELECT 1120 # undef zlup 1121 # undef zldw 1122 ! 1123 CASE ( 3 ) ! Bougeault & Lacarrere 89 length-scale 1124 ! 1125 # define zlup zRH 1126 # define zldw zFC 1127 ! zCF is used for matrix inversion 1128 ! 1129 DO jj = 1, jpj ! outer loop 1130 ! 1131 DO jk = 2, jpkam1 1132 DO ji = 1,jpi 1133 zcff = 1.0_wp / e3w_abl( jk )**2 1134 zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 1135 zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 1136 zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 1137 ENDDO 1138 ENDDO 1139 zsh2(:, 1) = zsh2( :, 2) 1140 zsh2(:, jpka) = zsh2( :, jpkam1) 1141 1142 DO ji = 1, jpi 1143 zcff = zrough(ji,jj) * rn_Lsfc 1144 zlup( ji, 1 ) = zcff 1145 zldw( ji, 1 ) = zcff 1146 zlup( ji, jpka ) = mxl_min 1147 zldw( ji, jpka ) = mxl_min 1148 END DO 1149 1150 DO jk = 2,jpka-1 1151 DO ji = 1, jpi 1152 zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 1153 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 1154 END DO 1155 END DO 1156 !! 1157 !! BL89 search for lup 1158 !! ---------------------------------------------------------- 1159 DO jk=2,jpka-1 1160 ! 1161 DO ji = 1, jpi 1162 zCF(ji,1:jpka) = 0._wp 1163 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 1164 ln_foundl(ji ) = .false. 1165 END DO 1166 ! 1167 DO jkup=jk+1,jpka-1 1168 DO ji = 1, jpi 1169 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1170 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 1171 zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 1172 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & 1173 & + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) & 1174 & + 0.5_wp * e3t_abl(jkup) * rn_Rod * & 1175 & ( SQRT(tke_abl( ji, jj, jkup , nt_a ))*zsh2(ji,jkup ) & 1176 & + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) ) 1177 1178 IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1179 zcff2 = ghw_abl(jkup ) - ghw_abl(jk) 1180 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 1181 zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / & 1182 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 1183 zlup(ji,jk) = zcff 1184 zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk) 1185 ln_foundl(ji) = .true. 1186 END IF 1187 END DO 1188 END DO 1189 ! 1190 END DO 1191 !! 1192 !! BL89 search for ldwn 1193 !! ---------------------------------------------------------- 1194 DO jk=2,jpka-1 1195 ! 1196 DO ji = 1, jpi 1197 zCF(ji,1:jpka) = 0._wp 1198 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 1199 ln_foundl(ji ) = .false. 1200 END DO 1201 ! 1202 DO jkdwn=jk-1,1,-1 1203 DO ji = 1, jpi 1204 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1205 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 1206 zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) & 1207 & * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 1208 & +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) & 1209 & + 0.5_wp * e3t_abl(jkup) * rn_Rod * & 1210 & ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) & 1211 & + SQRT(tke_abl( ji, jj, jkdwn , nt_a ))*zsh2(ji,jkdwn ) ) 1212 1213 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1214 zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 1215 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 1216 zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / & 1217 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1218 zldw(ji,jk) = zcff 1219 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1220 ln_foundl(ji) = .true. 1221 END IF 1222 END DO 1223 END DO 1224 ! 1225 END DO 1226 1227 DO jk = 1, jpka 1228 DO ji = 1, jpi 1229 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1230 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1231 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1232 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1233 END DO 1234 END DO 1235 1236 END DO 1237 # undef zlup 1238 # undef zldw 1239 ! 1240 ! 1241 END SELECT 934 1242 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 935 1243 ! ! Finalize the computation of turbulent visc./diff. 936 1244 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 937 1245 938 1246 !------------- 939 1247 DO jj = 1, jpj ! outer loop 940 1248 !------------- 941 DO jk = 1, jpka 1249 DO jk = 1, jpka 942 1250 DO ji = 1, jpi ! vector opt. 943 zcff = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk ) &944 & * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) )945 zcff2 946 zcff = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )947 !!FL: MAX function probably useless because of the definition of mxl_min 1251 zcff = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk ) & 1252 & * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) ) 1253 zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z) 1254 zcff = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) ) 1255 !!FL: MAX function probably useless because of the definition of mxl_min 948 1256 Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak ) 949 Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak ) 950 END DO 951 END DO 952 !------------- 953 END DO 1257 Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak ) 1258 END DO 1259 END DO 1260 !------------- 1261 END DO 954 1262 !------------- 955 1263 … … 969 1277 !! 970 1278 !! --------------------------------------------------------------------- 971 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk 972 1279 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk 1280 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d 973 1281 INTEGER :: ji,jj 974 975 976 1282 REAL(wp) :: smth_a, smth_b 1283 REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY 1284 REAL(wp) :: zumsk,zvmsk 977 1285 !! 978 1286 !!========================================================= … … 986 1294 zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk 987 1295 END_2D 988 989 1296 1297 DO_2D_10_11 990 1298 zvmsk = msk(ji,jj) * msk(ji,jj+1) 991 1299 zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk 992 993 994 1300 END_2D 1301 1302 DO_2D_10_00 995 1303 zFY ( ji, jj ) = zdY ( ji, jj ) & 996 1304 & + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) & 997 1305 & - (zdX ( ji, jj ) - zdX( ji-1, jj )) ) 998 1306 END_2D 999 1307 1000 1308 DO_2D_00_10 … … 1010 1318 & +zFY( ji, jj ) - zFY( ji, jj-1 ) ) 1011 1319 END_2D 1012 !! 1320 1013 1321 !--------------------------------------------------------------------------------------------------- 1014 1322 END SUBROUTINE smooth_pblh -
NEMO/trunk/src/ABL/ablrst.F90
r12649 r13214 109 109 CALL iom_delay_rst( 'WRITE', 'ABL', numraw ) ! save only abl delayed global communication variables 110 110 111 ! Prognostic variables111 ! Prognostic (after timestep + swap time indices = now timestep) variables 112 112 CALL iom_rstput( iter, nitrst, numraw, 'u_abl', u_abl(:,:,:,nt_n ) ) 113 113 CALL iom_rstput( iter, nitrst, numraw, 'v_abl', v_abl(:,:,:,nt_n ) ) … … 117 117 CALL iom_rstput( iter, nitrst, numraw, 'avm_abl', avm_abl(:,:,: ) ) 118 118 CALL iom_rstput( iter, nitrst, numraw, 'avt_abl', avt_abl(:,:,: ) ) 119 CALL iom_rstput( iter, nitrst, numraw, 'mxl_abl', mxl_abl(:,:,: ) )119 CALL iom_rstput( iter, nitrst, numraw,'mxld_abl',mxld_abl(:,:,: ) ) 120 120 CALL iom_rstput( iter, nitrst, numraw, 'pblh', pblh(:,: ) ) 121 121 ! … … 172 172 CALL iom_get( numrar, jpdom_autoglo, 'avm_abl', avm_abl(:,:,: ) ) 173 173 CALL iom_get( numrar, jpdom_autoglo, 'avt_abl', avt_abl(:,:,: ) ) 174 CALL iom_get( numrar, jpdom_autoglo, 'mxl_abl', mxl_abl(:,:,: ) )174 CALL iom_get( numrar, jpdom_autoglo,'mxld_abl',mxld_abl(:,:,: ) ) 175 175 CALL iom_get( numrar, jpdom_autoglo, 'pblh', pblh(:,: ) ) 176 176 CALL iom_delay_rst( 'READ', 'ABL', numrar ) ! read only abl delayed global communication variables -
NEMO/trunk/src/ABL/par_abl.F90
r12814 r13214 28 28 LOGICAL , PUBLIC :: ln_hpgls_frc !: forcing of ABL winds by large-scale pressure gradient 29 29 LOGICAL , PUBLIC :: ln_smth_pblh !: smoothing of atmospheric PBL height 30 !LOGICAL , PUBLIC :: ln_topbc_neumann = .FALSE. !: idealised testcases only 30 31 31 LOGICAL , PUBLIC :: ln_rstart_abl !: (de)activate abl restart32 CHARACTER(len=256), PUBLIC :: cn_ablrst_in !: suffix of abl restart name (input)33 CHARACTER(len=256), PUBLIC :: cn_ablrst_out !: suffix of abl restart name (output)34 CHARACTER(len=256), PUBLIC :: cn_ablrst_indir !: abl restart input directory35 CHARACTER(len=256), PUBLIC :: cn_ablrst_outdir !: abl restart output directory32 LOGICAL , PUBLIC :: ln_rstart_abl !: (de)activate abl restart 33 CHARACTER(len=256), PUBLIC :: cn_ablrst_in !: suffix of abl restart name (input) 34 CHARACTER(len=256), PUBLIC :: cn_ablrst_out !: suffix of abl restart name (output) 35 CHARACTER(len=256), PUBLIC :: cn_ablrst_indir !: abl restart input directory 36 CHARACTER(len=256), PUBLIC :: cn_ablrst_outdir !: abl restart output directory 36 37 37 38 !!--------------------------------------------------------------------- … … 46 47 REAL(wp), PUBLIC, PARAMETER :: rn_Cek = 258._wp !: Ekman constant for Richardson number 47 48 REAL(wp), PUBLIC, PARAMETER :: rn_epssfc = 1._wp / ( 1._wp + 2.8_wp * 2.8_wp ) 48 REAL(wp), PUBLIC :: rn_ ceps !: namelist parameter49 REAL(wp), PUBLIC :: rn_ cm !: namelist parameter50 REAL(wp), PUBLIC :: rn_ ct !: namelist parameter51 REAL(wp), PUBLIC :: rn_ ce !: namelist parameter49 REAL(wp), PUBLIC :: rn_Ceps !: namelist parameter 50 REAL(wp), PUBLIC :: rn_Cm !: namelist parameter 51 REAL(wp), PUBLIC :: rn_Ct !: namelist parameter 52 REAL(wp), PUBLIC :: rn_Ce !: namelist parameter 52 53 REAL(wp), PUBLIC :: rn_Rod !: namelist parameter 53 54 REAL(wp), PUBLIC :: rn_Sch 55 REAL(wp), PUBLIC :: rn_Esfc 56 REAL(wp), PUBLIC :: rn_Lsfc 54 57 REAL(wp), PUBLIC :: mxl_min 55 58 REAL(wp), PUBLIC :: rn_ldyn_min !: namelist parameter -
NEMO/trunk/src/ABL/sbcabl.F90
r13208 r13214 71 71 & ln_hpgls_frc, ln_geos_winds, nn_dyn_restore, & 72 72 & rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max, & 73 & nn_amxl, rn_ cm, rn_ct, rn_ce, rn_ceps, rn_Rod, rn_Ric, &73 & nn_amxl, rn_Cm, rn_Ct, rn_Ce, rn_Ceps, rn_Rod, rn_Ric, & 74 74 & ln_smth_pblh 75 75 !!--------------------------------------------------------------------- 76 76 77 ! Namelist namsbc_abl in reference namelist : ABL parameters77 ! Namelist namsbc_abl in reference namelist : ABL parameters 78 78 READ ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 ) 79 79 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist' ) 80 ! Namelist namsbc_abl in configuration namelist : ABL parameters 80 ! 81 ! Namelist namsbc_abl in configuration namelist : ABL parameters 81 82 READ ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 ) 82 83 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist' ) … … 165 166 rn_Sch = rn_ce / rn_cm 166 167 mxl_min = (avm_bak / rn_cm) / sqrt( tke_min ) 168 rn_Esfc = 1._wp / SQRT(rn_cm*rn_ceps) 169 rn_Lsfc = vkarmn * SQRT(SQRT(rn_cm*rn_ceps)) / rn_cm 167 170 168 171 IF(lwp) THEN … … 171 174 WRITE(numout,*) ' ~~~~~~~~~~~' 172 175 IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 173 IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height ' 176 IF(nn_amxl==1) WRITE(numout,*) 'Modified Deardorff 80 length-scale ' 177 IF(nn_amxl==2) WRITE(numout,*) 'Bougeault and Lacarrere length-scale ' 178 IF(nn_amxl==3) WRITE(numout,*) 'Rodier et al. length-scale ' 174 179 WRITE(numout,*) ' Minimum value of atmospheric TKE = ',tke_min,' m^2 s^-2' 175 180 WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m' … … 178 183 WRITE(numout,*) ' Constant for Schmidt number = ',rn_Sch 179 184 WRITE(numout,*) ' Constant for TKE dissipation = ',rn_Ceps 185 WRITE(numout,*) ' Constant for TKE sfc boundary condition = ',rn_Esfc 186 WRITE(numout,*) ' Constant for mxl sfc boundary condition = ',rn_Lsfc 180 187 END IF 181 188 … … 202 209 ! ABL timestep 203 210 rDt_abl = nn_fsbc * rn_Dt 211 IF(lwp) WRITE(numout,*) ' ABL timestep = ', rDt_abl,' s' 204 212 205 213 ! Check parameters for dynamics … … 248 256 zcff = 2._wp * omega * SIN( rad * 90._wp ) !++ fmax 249 257 rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 250 !!GS: alternative shape251 !rest_eq(:,:) = SIN( 0.5_wp*rpi*(zcff - ABS(ff_t(:,:))) / (zcff - 3.e-5) )**8252 !WHERE(ABS(ff_t(:,:)).LE.3.e-5) rest_eq(:,:) = 1._wp253 258 ELSE 254 259 rest_eq(:,:) = 1._wp … … 271 276 CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 272 277 273 u_abl(:,:,:,nt_n ) = sf(jp_wndi)%fnow(:,:,:)274 v_abl(:,:,:,nt_n ) = sf(jp_wndj)%fnow(:,:,:)278 u_abl(:,:,:,nt_n ) = sf(jp_wndi)%fnow(:,:,:) 279 v_abl(:,:,:,nt_n ) = sf(jp_wndj)%fnow(:,:,:) 275 280 tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 276 281 tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) … … 279 284 avm_abl(:,:,: ) = avm_bak 280 285 avt_abl(:,:,: ) = avt_bak 281 mxl_abl(:,:,: ) = mxl_min282 286 pblh (:,: ) = ghw_abl( 3 ) !<-- assume that the pbl contains 3 grid points 283 287 u_abl (:,:,:,nt_a ) = 0._wp … … 285 289 tq_abl (:,:,:,nt_a,: ) = 0._wp 286 290 tke_abl(:,:,:,nt_a ) = 0._wp 291 292 mxlm_abl(:,:,: ) = mxl_min 293 mxld_abl(:,:,: ) = mxl_min 287 294 ENDIF 288 295
Note: See TracChangeset
for help on using the changeset viewer.
