Changeset 12588
- Timestamp:
- 2020-03-23T18:21:59+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src
- Files:
-
- 7 edited
-
ABL/abl.F90 (modified) (3 diffs)
-
ABL/ablmod.F90 (modified) (35 diffs)
-
ABL/ablrst.F90 (modified) (2 diffs)
-
ABL/par_abl.F90 (modified) (2 diffs)
-
ABL/sbcabl.F90 (modified) (7 diffs)
-
ICE/icesbc.F90 (modified) (1 diff)
-
OCE/SBC/sbcblk.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/abl.F90
r12489 r12588 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] … … 38 39 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: msk_abl 39 40 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: rest_eq 41 42 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: cft_abl 43 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,: ) :: taux_abl 44 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,: ) :: tauy_abl 40 45 ! 41 46 INTEGER , PUBLIC :: nt_n, nt_a !: now / after indices (equal 1 or 2) … … 55 60 !!---------------------------------------------------------------------- 56 61 ! 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 ) 62 ALLOCATE( u_abl (1:jpi,1:jpj,1:jpka,jptime ), & 63 & v_abl (1:jpi,1:jpj,1:jpka,jptime ), & 64 & tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), & 65 & tke_abl (1:jpi,1:jpj,1:jpka,jptime ), & 66 & avm_abl (1:jpi,1:jpj,1:jpka ), & 67 & avt_abl (1:jpi,1:jpj,1:jpka ), & 68 & mxld_abl(1:jpi,1:jpj,1:jpka ), & 69 & mxlm_abl(1:jpi,1:jpj,1:jpka ), & 70 & cft_abl (1:jpi,1:jpj,1:jpka ), & 71 & fft_abl (1:jpi,1:jpj ), & 72 & pblh (1:jpi,1:jpj ), & 73 & taux_abl(1:jpi,1:jpj ), & 74 & tauy_abl(1:jpi,1:jpj ), & 75 & msk_abl (1:jpi,1:jpj ), & 76 & rest_eq (1:jpi,1:jpj ), & 77 & e3t_abl (1:jpka), e3w_abl(1:jpka) , & 78 & ght_abl (1:jpka), ghw_abl(1:jpka) , STAT=ierr ) 69 79 ! 70 80 abl_alloc = ierr -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90
r12489 r12588 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 20 USE sbcblk ! use rn_?fac … … 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" … … 39 39 !=================================================================================================== 40 40 SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, & ! in 41 & pu_dta, pv_dta, pt_dta, pq_dta, & 41 & pu_dta, pv_dta, pt_dta, pq_dta, & 42 42 & pslp_dta, pgu_dta, pgv_dta, & 43 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 & ) 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 101 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point) 102 #endif 103 103 ! 104 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 105 REAL(wp), DIMENSION(1:jpi,2:jpka ) :: zCF 107 106 ! 108 107 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a … … 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 = .FALSE. 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(Cdn_oce(ji,jj)) ) !<-- 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 175 DO ji = 2, jpim1 176 !DO ji = 1,jpi !!GS: to be checked 174 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 … … 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 * zztmp2 188 #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 209 DO ji = 1,jpi 207 210 zcff = 1._wp / z_elem_b( ji, 2 ) 208 211 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 209 212 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 213 END DO 214 215 DO jk = 3, jpka 216 DO ji = 1,jpi 214 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 ) … … 218 221 END DO 219 222 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 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 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 240 IF( SemiImp_Cor ) THEN 241 242 !------------- 239 243 DO jk = 2, jpka ! outer loop 240 244 !------------- 241 ! 245 ! 242 246 ! Advance u_abl & v_abl to time n+1 243 247 DO_2D_11_11 244 248 zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl ) ! (f dt)**2 245 249 246 250 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 247 251 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n ) & 248 252 & + rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj , jk, nt_n ) ) & 249 253 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 250 254 251 255 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 252 256 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n ) & 253 257 & - rDt_abl * fft_abl(ji, jj) * u_abl ( ji , jj, jk, nt_n ) ) & 254 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 258 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 255 259 END_2D 256 ! 260 ! 257 261 !------------- 258 262 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 !------------- 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 ELSE 279 280 !------------- 281 DO jk = 2, jpka ! outer loop 282 !------------- 283 ! 284 IF( MOD( kt, 2 ) == 0 ) then 285 ! Advance u_abl & v_abl to time n+1 286 DO jj = 1, jpj 287 DO ji = 1, jpi 288 zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_n ) - pgv_dta(ji ,jj ,jk) ) 289 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff 290 zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_a ) - pgu_dta(ji ,jj ,jk) ) 291 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff ) 292 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * u_abl( ji, jj, jk, nt_a ) 293 END DO 294 END DO 295 ! 296 ELSE 297 ! 298 DO jj = 1, jpj 299 DO ji = 1, jpi 300 zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj , jk, nt_n ) - pgu_dta(ji ,jj ,jk) ) 301 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff 302 zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj , jk, nt_a ) - pgv_dta(ji ,jj ,jk) ) 303 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff ) 304 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) * v_abl( ji, jj, jk, nt_a ) 305 END DO 306 END DO 307 ! 308 ENDIF 309 ! 310 !------------- 311 END DO ! end outer loop !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 312 !------------- 313 314 ENDIF 315 316 !------------- 317 ! 318 IF( ln_hpgls_frc ) THEN 319 DO jj = 1, jpj ! outer loop 263 320 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 321 DO ji = 1, jpi 279 322 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) 323 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 281 324 ENDDO 282 325 ENDDO 283 ENDDO 326 ENDDO 284 327 END IF 285 286 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 287 ! ! 4 *** Advance u,v to time n+1 288 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 289 ! 290 ! Vertical diffusion for u_abl 328 329 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 330 ! ! 4 *** Advance u,v to time n+1 331 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 332 ! 333 ! Vertical diffusion for u_abl 291 334 !------------- 292 335 DO jj = 1, jpj ! outer loop 293 !------------- 336 !------------- 294 337 295 338 DO jk = 3, jpkam1 296 DO ji = 1, jpi 339 DO ji = 1, jpi 297 340 z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 298 z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 341 z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 299 342 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 300 343 END DO 301 END DO 302 303 DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi) 344 END DO 345 346 DO ji = 2, jpi ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi) 304 347 !++ Surface boundary condition 305 348 z_elem_a( ji, 2 ) = 0._wp 306 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 349 z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 350 ! 351 zztmp1 = pcd_du(ji, jj) 352 zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) * rn_vfac 353 #if defined key_si3 311 354 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 355 zzice = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) ) * rn_vfac 356 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 357 #endif 358 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 316 359 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 360 361 IF( ln_topbc_neumann ) THEN 362 !++ Top Neumann B.C. 363 z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 364 z_elem_c( ji, jpka ) = 0._wp 365 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 366 !u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl ( ji, jj, jpka, nt_a ) 367 ELSE 368 !++ Top Dirichlet B.C. 369 z_elem_a( ji, jpka ) = 0._wp 370 z_elem_c( ji, jpka ) = 0._wp 371 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 372 u_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 373 ENDIF 374 375 END DO 328 376 !! 329 377 !! Matrix inversion 330 378 !! ---------------------------------------------------------- 331 DO ji = 2, jpi 379 DO ji = 2, jpi 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 411 DO ji = 1, jpi 364 412 z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 365 z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 413 z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 366 414 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 367 415 END DO 368 416 END DO 369 417 370 DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 418 DO ji = 1, jpi ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 371 419 !++ Surface boundary condition 372 420 z_elem_a( ji, 2 ) = 0._wp 373 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 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 IF( ln_topbc_neumann ) THEN 434 !++ Top Neumann B.C. 435 z_elem_a( ji, jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 436 z_elem_c( ji, jpka ) = 0._wp 437 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 438 !v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl ( ji, jj, jpka, nt_a ) 439 ELSE 440 !++ Top Dirichlet B.C. 441 z_elem_a( ji, jpka ) = 0._wp 442 z_elem_c( ji, jpka ) = 0._wp 443 z_elem_b( ji, jpka ) = e3t_abl( jpka ) 444 v_abl ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 445 ENDIF 446 447 END DO 394 448 !! 395 449 !! Matrix inversion 396 450 !! ---------------------------------------------------------- 397 DO ji = 1, jpi 451 DO ji = 1, jpi 398 452 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 453 zCF (ji, 2 ) = - zcff * z_elem_c( ji, 2 ) 454 v_abl (ji,jj,2,nt_a) = zcff * v_abl ( ji, jj, 2, nt_a ) 455 END DO 456 457 DO jk = 3, jpka 458 DO ji = 1, jpi 405 459 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF (ji, jk-1 ) ) 406 460 zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) … … 409 463 END DO 410 464 END DO 411 412 DO jk = jpkam1,2,-1 413 DO ji = 1, jpi 465 466 DO jk = jpkam1,2,-1 467 DO ji = 1, jpi 414 468 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 469 END DO 416 470 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 471 ! 472 !------------- 473 END DO ! end outer loop 474 !------------- 475 476 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 477 ! ! 5 *** Apply nudging on the dynamics and the tracers 478 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 479 427 480 IF( nn_dyn_restore > 0 ) THEN 428 !------------- 481 !------------- 429 482 DO jk = 2, jpka ! outer loop 430 !------------- 483 !------------- 431 484 DO_2D_01_01 432 485 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) ) 486 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 487 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 435 488 zmsk = msk_abl(ji,jj) 436 489 zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2 & 437 490 & + jp_alp1_dyn * zsig + jp_alp0_dyn 438 491 zcff = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt ! zcff = 1 for masked points 439 ! rn_Dt = rDt_abl / nn_fsbc 492 ! rn_Dt = rDt_abl / nn_fsbc 440 493 zcff = zcff * rest_eq(ji,jj) 441 z_cft( ji, jj, jk ) = zcff442 494 u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) & 443 & + zcff * pu_dta( ji, jj, jk ) 495 & + zcff * pu_dta( ji, jj, jk ) 444 496 v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * v_abl( ji, jj, jk, nt_a ) & 445 497 & + zcff * pv_dta( ji, jj, jk ) … … 447 499 !------------- 448 500 END DO ! end outer loop 449 !------------- 501 !------------- 450 502 END IF 451 503 452 !------------- 504 !------------- 453 505 DO jk = 2, jpka ! outer loop 454 !------------- 506 !------------- 455 507 DO_2D_11_11 456 508 zcff1 = pblh( ji, jj ) 457 509 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 458 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 510 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 459 511 zmsk = msk_abl(ji,jj) 460 512 zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2 & 461 513 & + jp_alp1_tra * zsig + jp_alp0_tra 462 514 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 515 ! rn_Dt = rDt_abl / nn_fsbc 465 516 tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) & 466 517 & + zcff * pt_dta( ji, jj, jk ) 467 518 468 519 tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa ) & 469 520 & + zcff * pq_dta( ji, jj, jk ) 470 521 471 522 END_2D 472 523 !------------- 473 524 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.)525 !------------- 526 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 527 ! ! 6 *** MPI exchanges 528 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 529 ! 530 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1., v_abl(:,:,:,nt_a) , 'T', -1., kfillmode = jpfillnothing ) 480 531 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 532 ! 533 #if defined key_iomput 482 534 ! first ABL level 483 535 IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl", u_abl(:,:,2,nt_a ) ) … … 497 549 IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka ) ) 498 550 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) )551 IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxlm_abl(:,:,2:jpka ) ) 500 552 IF ( iom_use("pblh" ) ) CALL iom_put ( "pblh" , pblh(:,: ) ) 501 553 ! debug (to be removed) … … 504 556 IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta", pt_dta(:,:,2:jpka) ) 505 557 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 558 IF( ln_geos_winds ) THEN 508 559 IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2 ) ) … … 513 564 IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 514 565 END IF 515 ! 566 #endif 516 567 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 517 568 ! ! 7 *** Finalize flux computation 518 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 519 569 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 570 ! 520 571 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 572 ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta ) 573 zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa ) 574 zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 575 psen( ji, jj ) = - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) !GS: negative sign to respect aerobulk convention 576 pevp( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) ) 577 rhoa( ji, jj ) = zcff 529 578 END_2D 530 579 531 580 DO_2D_01_01 532 zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji ,jj) + pssu(ji-1,jj) )533 zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj ) + pssv(ji,jj-1) )581 zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji,jj) + pssu(ji-1,jj ) ) 582 zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj) + pssv(ji ,jj-1) ) 534 583 END_2D 535 ! 584 ! 536 585 CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. ) 537 586 ! 538 587 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 539 588 DO_2D_11_11 540 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 541 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 542 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 543 544 pwndm (ji,jj) = zcff 545 ptaum (ji,jj) = zztmp * zcff 546 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 547 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 589 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 590 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 591 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 592 pwndm (ji,jj) = zcff 593 ptaum (ji,jj) = zztmp * zcff 594 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 595 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 596 taux_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_i(ji,jj) 597 tauy_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_j(ji,jj) 548 598 END_2D 549 599 ! ... utau, vtau at U- and V_points, resp. … … 574 624 ! ------------------------------------------------------------ ! 575 625 DO_2D_00_00 576 577 zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) )578 zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )579 626 627 zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 628 zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 629 580 630 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 581 631 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & … … 592 642 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 593 643 ! ! 8 *** Swap time indices for the next timestep 594 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 644 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 595 645 nt_n = 1 + MOD( kt , 2) 596 646 nt_a = 1 + MOD( kt+1, 2) 597 ! 647 ! 598 648 !--------------------------------------------------------------------------------------------------- 599 649 END SUBROUTINE abl_stp … … 634 684 !! (= Kz dz[Ub] * dz[Un] ) 635 685 !! --------------------------------------------------------------------- 636 INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn 686 INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn 637 687 INTEGER, DIMENSION(1:jpi ) :: ikbl 638 688 REAL(wp) :: zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv 639 REAL(wp) :: zdU , zdV, zcff1,zshear,zbuoy,zsig, zustar2640 REAL(wp) :: zdU2, zdV2641 REAL(wp) :: zwndi, zwndj689 REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2 690 REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 691 REAL(wp) :: zwndi, zwndj 642 692 REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2 643 693 REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2 644 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF 694 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF 645 695 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a 646 696 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b … … 648 698 LOGICAL :: ln_Patankar = .FALSE. 649 699 LOGICAL :: ln_dumpvar = .FALSE. 650 LOGICAL , DIMENSION(1:jpi ) :: ln_foundl 700 LOGICAL , DIMENSION(1:jpi ) :: ln_foundl 651 701 ! 652 702 tind = nt_n … … 660 710 !------------- 661 711 ! 662 ! Compute vertical shear 712 ! Compute vertical shear 663 713 DO jk = 2, jpkam1 664 DO ji = 1,jpi 665 zcff = 1.0_wp / e3w_abl( jk )**2 666 zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 714 DO ji = 1,jpi 715 zcff = 1.0_wp / e3w_abl( jk )**2 716 zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 667 717 zdV = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 668 zsh2(ji,jk) = zdU+zdV 669 END DO 670 END DO 718 zsh2(ji,jk) = zdU+zdV !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 ) 719 END DO 720 END DO 671 721 ! 672 722 ! Compute brunt-vaisala frequency 673 723 DO jk = 2, jpkam1 674 DO ji = 1,jpi 675 zcff = grav * itvref / e3w_abl( jk ) 724 DO ji = 1,jpi 725 zcff = grav * itvref / e3w_abl( jk ) 676 726 zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta) 677 727 zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) & … … 679 729 zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi) 680 730 END DO 681 END DO 731 END DO 682 732 ! 683 733 ! Terms for the tridiagonal problem 684 734 DO jk = 2, jpkam1 685 DO ji = 1,jpi 735 DO ji = 1,jpi 686 736 zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point 687 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 688 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 689 737 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 738 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 739 690 740 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 691 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 741 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 692 742 IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE 693 743 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 694 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl _abl(ji,jj,jk) ! diagonal695 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )! right-hand-side744 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) ! diagonal 745 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) ) ! right-hand-side 696 746 ELSE 697 747 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 698 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl _abl(ji,jj,jk) & ! diagonal699 & - e3w_abl(jk) * rDt_abl * zbuoy 700 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side 748 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) & ! diagonal 749 & - e3w_abl(jk) * rDt_abl * zbuoy 750 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side 701 751 END IF 702 752 END DO 703 END DO 704 705 DO ji = 1,jpi ! vector opt. 706 zesrf = MAX( 4.63_wp * ustar2(ji,jj), tke_min ) 707 zetop = tke_min 708 z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp 709 z_elem_a ( ji, jpka ) = 0._wp; z_elem_c ( ji, jpka ) = 0._wp; z_elem_b ( ji, jpka ) = 1._wp 710 tke_abl( ji, jj, 1, nt_a ) = zesrf 711 tke_abl( ji, jj, jpka, nt_a ) = zetop 712 zbn2(ji,jj, 1) = zbn2( ji,jj, 2) 713 zsh2(ji, 1) = zsh2( ji, 2) 714 zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 753 END DO 754 755 DO ji = 1,jpi ! vector opt. 756 zesrf = MAX( rn_esfc * ustar2(ji,jj), tke_min ) 757 zetop = tke_min 758 759 z_elem_a ( ji, 1 ) = 0._wp; z_elem_c ( ji, 1 ) = 0._wp; z_elem_b ( ji, 1 ) = 1._wp 760 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 ) 761 z_elem_c ( ji, jpka ) = 0._wp 762 z_elem_b ( ji, jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 763 764 tke_abl ( ji, jj, 1 , nt_a ) = zesrf 765 tke_abl ( ji, jj, jpka , nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 766 767 zbn2(ji,jj, 1) = zbn2( ji,jj, 2) 768 zsh2(ji, 1) = zsh2( ji, 2) 769 zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 715 770 zsh2(ji, jpka) = zsh2( ji , jpkam1) 716 END DO 771 END DO 717 772 !! 718 773 !! Matrix inversion … … 720 775 DO ji = 1,jpi 721 776 zcff = 1._wp / z_elem_b( ji, 1 ) 722 zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 ) 723 tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a ) 724 END DO 725 726 DO jk = 2, jpka 777 zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 ) 778 tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a ) 779 END DO 780 781 DO jk = 2, jpka 727 782 DO ji = 1,jpi 728 783 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) ) … … 732 787 END DO 733 788 END DO 734 735 DO jk = jpkam1,1,-1 789 790 DO jk = jpkam1,1,-1 736 791 DO ji = 1,jpi 737 792 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) 738 793 END DO 739 794 END DO 740 741 !!FL should not be needed because of Patankar procedure 795 796 !!FL should not be needed because of Patankar procedure 742 797 tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min ) 743 798 … … 745 800 !! Diagnose PBL height 746 801 !! ---------------------------------------------------------- 747 748 749 ! 802 803 804 ! 750 805 ! arrays zRH, zFC and zCF are available at this point 751 806 ! and zFC(:, 1 ) = 0. 752 807 ! diagnose PBL height based on zsh2 and zbn2 753 808 zFC ( : ,1) = 0._wp 754 ikbl( 1:jpi ) = 0 755 809 ikbl( 1:jpi ) = 0 810 756 811 DO jk = 2,jpka 757 DO ji = 1, jpi 812 DO ji = 1, jpi 758 813 zcff = ghw_abl( jk-1 ) 759 814 zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) ) … … 781 836 ELSE 782 837 pblh( ji, jj ) = ghw_abl(jpka) 783 END IF 784 END DO 785 !------------- 786 END DO 787 !------------- 788 ! 789 ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 838 END IF 839 END DO 840 !------------- 841 END DO 842 !------------- 843 ! 844 ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 790 845 IF(ln_smth_pblh) THEN 791 CALL lbc_lnk( 'ablmod', pblh, 'T', 1. )846 CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing) 792 847 CALL smooth_pblh( pblh, msk_abl ) 793 CALL lbc_lnk( 'ablmod', pblh, 'T', 1. )848 CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing) 794 849 ENDIF 795 850 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 799 854 SELECT CASE ( nn_amxl ) 800 855 ! 801 CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom 802 # define zlup zRH 803 # define zldw zFC 856 CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom 857 # define zlup zRH 858 # define zldw zFC 804 859 DO jj = 1, jpj ! outer loop 805 860 ! 806 861 DO ji = 1, jpi 807 mxl_abl ( ji, jj, 1 ) = 0._wp 808 mxl_abl ( ji, jj, jpka ) = mxl_min 809 zldw( ji, 1 ) = 0._wp 810 zlup( ji, jpka ) = 0._wp 811 END DO 812 ! 813 DO jk = 2, jpkam1 814 DO ji = 1, jpi 815 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 816 mxl_abl( ji, jj, jk ) = MAX( mxl_min, & 817 & SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) ) 818 END DO 819 END DO 862 mxld_abl ( ji, jj, 1 ) = 0._wp 863 mxld_abl ( ji, jj, jpka ) = mxl_min 864 mxlm_abl ( ji, jj, 1 ) = 0._wp 865 mxlm_abl ( ji, jj, jpka ) = mxl_min 866 zldw ( ji, 1 ) = zrough(ji,jj) * rn_Lsfc 867 zlup ( ji, jpka ) = mxl_min 868 END DO 869 ! 870 DO jk = 2, jpkam1 871 DO ji = 1, jpi 872 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 873 mxlm_abl( ji, jj, jk ) = MAX( mxl_min, & 874 & SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) ) 875 END DO 876 END DO 820 877 ! 821 878 ! Limit mxl 822 DO jk = jpkam1,1,-1 823 DO ji = 1, jpi 824 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl _abl(ji, jj, jk) )825 END DO 826 END DO 879 DO jk = jpkam1,1,-1 880 DO ji = 1, jpi 881 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 882 END DO 883 END DO 827 884 ! 828 885 DO jk = 2, jpka 829 DO ji = 1, jpi 830 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) ) 831 END DO 832 END DO 886 DO ji = 1, jpi 887 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 888 END DO 889 END DO 890 ! 891 ! DO jk = 1, jpka 892 ! DO ji = 1, jpi 893 ! mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 894 ! mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ), zlup( ji, jk ) ) 895 ! END DO 896 ! END DO 833 897 ! 834 898 DO jk = 1, jpka 835 899 DO ji = 1, jpi 836 mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 837 END DO 838 END DO 839 ! 840 END DO 841 # undef zlup 842 # undef zldw 843 ! 844 ! 845 CASE ( 1 ) ! length-scale computed as the distance to the PBL height 846 DO jj = 1,jpj ! outer loop 847 ! 848 DO ji = 1, jpi ! vector opt. 849 zcff = 1._wp / pblh( ji, jj ) ! inverse of hbl 850 DO jk = 1, jpka 851 zsig = MIN( zcff * ghw_abl( jk ), 1. ) 852 zcff1 = pblh( ji, jj ) 853 mxl_abl( ji, jj, jk ) = mxl_min & 854 & + zsig * ( amx1*zcff1 + bmx1*mxl_min ) & 855 & + zsig * zsig * ( amx2*zcff1 + bmx2*mxl_min ) & 856 & + zsig**3 * ( amx3*zcff1 + bmx3*mxl_min ) & 857 & + zsig**4 * ( amx4*zcff1 + bmx4*mxl_min ) & 858 & + zsig**5 * ( amx5*zcff1 + bmx5*mxl_min ) 859 END DO 860 END DO 861 ! 862 END DO 863 ! 900 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 901 ! zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 902 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 903 !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) ) 904 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 905 END DO 906 END DO 907 ! 908 END DO 909 # undef zlup 910 # undef zldw 911 ! 912 ! 913 CASE ( 1 ) ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom 914 # define zlup zRH 915 # define zldw zFC 916 DO jj = 1, jpj ! outer loop 917 ! 918 DO jk = 2, jpkam1 919 DO ji = 1,jpi 920 zcff = 1.0_wp / e3w_abl( jk )**2 921 zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 922 zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 923 zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 924 ENDDO 925 ENDDO 926 ! 927 DO ji = 1, jpi 928 zcff = zrough(ji,jj) * rn_Lsfc 929 mxld_abl ( ji, jj, 1 ) = zcff 930 mxld_abl ( ji, jj, jpka ) = mxl_min 931 mxlm_abl ( ji, jj, 1 ) = zcff 932 mxlm_abl ( ji, jj, jpka ) = mxl_min 933 zldw ( ji, 1 ) = zcff 934 zlup ( ji, jpka ) = mxl_min 935 END DO 936 ! 937 DO jk = 2, jpkam1 938 DO ji = 1, jpi 939 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 940 zcff = 2.*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) & 941 & + SQRT( rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.*zbuoy ) ) 942 mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff ) 943 END DO 944 END DO 945 ! 946 ! Limit mxl 947 DO jk = jpkam1,1,-1 948 DO ji = 1, jpi 949 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 950 END DO 951 END DO 952 ! 953 DO jk = 2, jpka 954 DO ji = 1, jpi 955 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 956 END DO 957 END DO 958 ! 959 DO jk = 1, jpka 960 DO ji = 1, jpi 961 !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 962 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 963 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 964 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 965 !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) ) 966 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 967 END DO 968 END DO 969 ! 970 END DO 971 # undef zlup 972 # undef zldw 973 ! ! 864 974 CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale 865 975 ! 866 # define zlup zRH 867 # define zldw zFC 976 # define zlup zRH 977 # define zldw zFC 868 978 ! zCF is used for matrix inversion 869 ! 979 ! 870 980 DO jj = 1, jpj ! outer loop 871 872 DO ji = 1, jpi 873 zlup( ji, 1 ) = mxl_min 874 zldw( ji, 1 ) = mxl_min 981 982 DO ji = 1, jpi 983 zcff = zrough(ji,jj) * rn_Lsfc 984 zlup( ji, 1 ) = zcff 985 zldw( ji, 1 ) = zcff 875 986 zlup( ji, jpka ) = mxl_min 876 zldw( ji, jpka ) = mxl_min 877 END DO 878 987 zldw( ji, jpka ) = mxl_min 988 END DO 989 879 990 DO jk = 2,jpka-1 880 991 DO ji = 1, jpi 881 992 zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 882 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 883 END DO 884 END DO 993 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 994 END DO 995 END DO 885 996 !! 886 997 !! BL89 search for lup 887 !! ---------------------------------------------------------- 888 DO jk=2,jpka-1 998 !! ---------------------------------------------------------- 999 DO jk=2,jpka-1 889 1000 ! 890 1001 DO ji = 1, jpi … … 892 1003 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 893 1004 ln_foundl(ji ) = .false. 894 END DO 895 ! 1005 END DO 1006 ! 896 1007 DO jkup=jk+1,jpka-1 897 1008 DO ji = 1, jpi 1009 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1010 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 898 1011 zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 899 & ( zb n2(ji,jj,jkup )*(ghw_abl(jkup )-ghw_abl(jk)) &900 & + zb n2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) )1012 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & 1013 & + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) 901 1014 IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 902 1015 zcff2 = ghw_abl(jkup ) - ghw_abl(jk) 903 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 1016 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 904 1017 zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / & 905 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 906 zlup(ji,jk) = zcff 1018 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 1019 zlup(ji,jk) = zcff 1020 zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk) 907 1021 ln_foundl(ji) = .true. 908 1022 END IF … … 910 1024 END DO 911 1025 ! 912 END DO 1026 END DO 913 1027 !! 914 1028 !! BL89 search for ldwn 915 !! ---------------------------------------------------------- 916 DO jk=2,jpka-1 1029 !! ---------------------------------------------------------- 1030 DO jk=2,jpka-1 917 1031 ! 918 1032 DO ji = 1, jpi … … 920 1034 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 921 1035 ln_foundl(ji ) = .false. 922 END DO 923 ! 1036 END DO 1037 ! 924 1038 DO jkdwn=jk-1,1,-1 925 DO ji = 1, jpi 1039 DO ji = 1, jpi 1040 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1041 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 926 1042 zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) & 927 & * ( zb n2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &928 + zb n2(ji,jj,jkdwn )*(ghw_abl(jk)-ghw_abl(jkdwn )) )929 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1043 & * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 1044 + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) 1045 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 930 1046 zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 931 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 1047 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 932 1048 zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / & 933 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 934 zldw(ji,jk) = zcff 935 ln_foundl(ji) = .true. 936 END IF 937 END DO 938 END DO 939 ! 1049 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1050 zldw(ji,jk) = zcff 1051 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1052 ln_foundl(ji) = .true. 1053 END IF 1054 END DO 1055 END DO 1056 ! 940 1057 END DO 941 1058 942 1059 DO jk = 1, jpka 943 DO ji = 1, jpi 944 mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min ) 945 END DO 946 END DO 1060 DO ji = 1, jpi 1061 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1062 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1063 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1064 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1065 END DO 1066 END DO 947 1067 948 1068 END DO 949 # undef zlup 950 # undef zldw 951 ! 952 END SELECT 1069 # undef zlup 1070 # undef zldw 1071 ! 1072 CASE ( 3 ) ! Bougeault & Lacarrere 89 length-scale 1073 ! 1074 # define zlup zRH 1075 # define zldw zFC 1076 ! zCF is used for matrix inversion 1077 ! 1078 DO jj = 1, jpj ! outer loop 1079 ! 1080 DO jk = 2, jpkam1 1081 DO ji = 1,jpi 1082 zcff = 1.0_wp / e3w_abl( jk )**2 1083 zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 1084 zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 1085 zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 1086 ENDDO 1087 ENDDO 1088 zsh2(:, 1) = zsh2( :, 2) 1089 zsh2(:, jpka) = zsh2( :, jpkam1) 1090 1091 DO ji = 1, jpi 1092 zcff = zrough(ji,jj) * rn_Lsfc 1093 zlup( ji, 1 ) = zcff 1094 zldw( ji, 1 ) = zcff 1095 zlup( ji, jpka ) = mxl_min 1096 zldw( ji, jpka ) = mxl_min 1097 END DO 1098 1099 DO jk = 2,jpka-1 1100 DO ji = 1, jpi 1101 zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 1102 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 1103 END DO 1104 END DO 1105 !! 1106 !! BL89 search for lup 1107 !! ---------------------------------------------------------- 1108 DO jk=2,jpka-1 1109 ! 1110 DO ji = 1, jpi 1111 zCF(ji,1:jpka) = 0._wp 1112 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 1113 ln_foundl(ji ) = .false. 1114 END DO 1115 ! 1116 DO jkup=jk+1,jpka-1 1117 DO ji = 1, jpi 1118 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1119 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 1120 zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 1121 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & 1122 & + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) & 1123 & + 0.5_wp * e3t_abl(jkup) * rn_Rod * & 1124 & ( SQRT(tke_abl( ji, jj, jkup , nt_a ))*zsh2(ji,jkup ) & 1125 & + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) ) 1126 1127 IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1128 zcff2 = ghw_abl(jkup ) - ghw_abl(jk) 1129 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 1130 zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / & 1131 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 1132 zlup(ji,jk) = zcff 1133 zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk) 1134 ln_foundl(ji) = .true. 1135 END IF 1136 END DO 1137 END DO 1138 ! 1139 END DO 1140 !! 1141 !! BL89 search for ldwn 1142 !! ---------------------------------------------------------- 1143 DO jk=2,jpka-1 1144 ! 1145 DO ji = 1, jpi 1146 zCF(ji,1:jpka) = 0._wp 1147 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 1148 ln_foundl(ji ) = .false. 1149 END DO 1150 ! 1151 DO jkdwn=jk-1,1,-1 1152 DO ji = 1, jpi 1153 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1154 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 1155 zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) & 1156 & * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 1157 & +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) & 1158 & + 0.5_wp * e3t_abl(jkup) * rn_Rod * & 1159 & ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) & 1160 & + SQRT(tke_abl( ji, jj, jkdwn , nt_a ))*zsh2(ji,jkdwn ) ) 1161 1162 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1163 zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 1164 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 1165 zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / & 1166 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1167 zldw(ji,jk) = zcff 1168 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1169 ln_foundl(ji) = .true. 1170 END IF 1171 END DO 1172 END DO 1173 ! 1174 END DO 1175 1176 DO jk = 1, jpka 1177 DO ji = 1, jpi 1178 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1179 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1180 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1181 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1182 END DO 1183 END DO 1184 1185 END DO 1186 # undef zlup 1187 # undef zldw 1188 ! 1189 ! 1190 ! 1191 END SELECT 953 1192 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 954 1193 ! ! Finalize the computation of turbulent visc./diff. 955 1194 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 956 1195 957 1196 !------------- 958 1197 DO jj = 1, jpj ! outer loop 959 1198 !------------- 960 DO jk = 1, jpka 1199 DO jk = 1, jpka 961 1200 DO ji = 1, jpi ! vector opt. 962 zcff = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk ) &963 & * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) )964 zcff2 965 zcff = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )966 !!FL: MAX function probably useless because of the definition of mxl_min 1201 zcff = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk ) & 1202 & * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) ) 1203 zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z) 1204 zcff = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) ) 1205 !!FL: MAX function probably useless because of the definition of mxl_min 967 1206 Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff , avm_bak ) 968 Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak ) 969 END DO 970 END DO 971 !------------- 972 END DO 1207 Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak ) 1208 END DO 1209 END DO 1210 !------------- 1211 END DO 973 1212 !------------- 974 1213 … … 988 1227 !! 989 1228 !! --------------------------------------------------------------------- 990 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk 991 1229 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk 1230 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d 992 1231 INTEGER :: ji,jj 993 994 995 1232 REAL(wp) :: smth_a, smth_b 1233 REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY 1234 REAL(wp) :: zumsk,zvmsk 996 1235 !! 997 1236 !!========================================================= … … 1005 1244 zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk 1006 1245 END_2D 1007 1008 1246 1247 DO_2D_10_11 1009 1248 zvmsk = msk(ji,jj) * msk(ji,jj+1) 1010 1249 zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk 1011 1012 1013 1250 END_2D 1251 1252 DO_2D_10_00 1014 1253 zFY ( ji, jj ) = zdY ( ji, jj ) & 1015 1254 & + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) & 1016 1255 & - (zdX ( ji, jj ) - zdX( ji-1, jj )) ) 1017 1256 END_2D 1018 1257 1019 1258 DO_2D_00_10 … … 1029 1268 & +zFY( ji, jj ) - zFY( ji, jj-1 ) ) 1030 1269 END_2D 1031 !! 1270 1032 1271 !--------------------------------------------------------------------------------------------------- 1033 1272 END SUBROUTINE smooth_pblh -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablrst.F90
r11945 r12588 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(:,:,: ) )120 119 CALL iom_rstput( iter, nitrst, numraw, 'pblh', pblh(:,: ) ) 121 120 ! … … 172 171 CALL iom_get( numrar, jpdom_autoglo, 'avm_abl', avm_abl(:,:,: ) ) 173 172 CALL iom_get( numrar, jpdom_autoglo, 'avt_abl', avt_abl(:,:,: ) ) 174 CALL iom_get( numrar, jpdom_autoglo, 'mxl_abl', mxl_abl(:,:,: ) )175 173 CALL iom_get( numrar, jpdom_autoglo, 'pblh', pblh(:,: ) ) 176 174 CALL iom_delay_rst( 'READ', 'ABL', numrar ) ! read only abl delayed global communication variables -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/par_abl.F90
r12489 r12588 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. !: smoothing of atmospheric PBL height 30 31 31 CHARACTER(len=256), PUBLIC :: cn_ablrst_in !: suffix of abl restart name (input) 32 CHARACTER(len=256), PUBLIC :: cn_ablrst_out !: suffix of abl restart name (output) 33 CHARACTER(len=256), PUBLIC :: cn_ablrst_indir !: abl restart input directory 34 CHARACTER(len=256), PUBLIC :: cn_ablrst_outdir !: abl restart output directory 32 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 35 37 36 38 !!--------------------------------------------------------------------- … … 45 47 REAL(wp), PUBLIC, PARAMETER :: rn_Cek = 258._wp !: Ekman constant for Richardson number 46 48 REAL(wp), PUBLIC, PARAMETER :: rn_epssfc = 1._wp / ( 1._wp + 2.8_wp * 2.8_wp ) 47 REAL(wp), PUBLIC :: rn_ ceps !: namelist parameter48 REAL(wp), PUBLIC :: rn_ cm !: namelist parameter49 REAL(wp), PUBLIC :: rn_ ct !: namelist parameter50 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 51 53 REAL(wp), PUBLIC :: rn_Rod !: namelist parameter 52 54 REAL(wp), PUBLIC :: rn_Sch 55 REAL(wp), PUBLIC :: rn_Esfc 56 REAL(wp), PUBLIC :: rn_Lsfc 53 57 REAL(wp), PUBLIC :: mxl_min 54 58 REAL(wp), PUBLIC :: rn_ldyn_min !: namelist parameter -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90
r12549 r12588 68 68 LOGICAL :: lluldl 69 69 NAMELIST/namsbc_abl/ cn_dir, cn_dom, cn_ablrst_in, cn_ablrst_out, & 70 & cn_ablrst_indir, cn_ablrst_outdir, 70 & cn_ablrst_indir, cn_ablrst_outdir, ln_rstart_abl, & 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 ) 167 168 rn_Esfc = 1._wp / SQRT(rn_cm*rn_ceps) 169 rn_Lsfc = vkarmn * SQRT(SQRT(rn_cm*rn_ceps)) / rn_cm 170 168 171 IF(lwp) THEN 169 172 WRITE(numout,*) … … 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==1) 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 … … 248 255 zcff = 2._wp * omega * SIN( rad * 90._wp ) !++ fmax 249 256 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 257 ELSE 254 258 rest_eq(:,:) = 1._wp … … 267 271 268 272 ! initialize ABL from data or restart 269 IF( ln_rstart ) THEN273 IF( ln_rstart_abl ) THEN 270 274 CALL abl_rst_read 271 275 ELSE 272 276 CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 273 277 274 u_abl(:,:,:,nt_n ) = sf(jp_wndi)%fnow(:,:,:)275 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(:,:,:) 276 280 tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 277 281 tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) … … 280 284 avm_abl(:,:,: ) = avm_bak 281 285 avt_abl(:,:,: ) = avt_bak 282 mxl_abl(:,:,: ) = mxl_min283 286 pblh (:,: ) = ghw_abl( 3 ) !<-- assume that the pbl contains 3 grid points 284 287 u_abl (:,:,:,nt_a ) = 0._wp -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icesbc.F90
r12377 r12588 148 148 CASE( jp_blk, jp_abl ) !--- bulk formulation & ABL formulation 149 149 CALL blk_ice_2 ( t_su, h_s, h_i, alb_ice, sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & 150 & sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) ) !150 & sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) ) 151 151 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 152 152 IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist ) -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90
r12489 r12588 104 104 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 105 105 ! 106 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice107 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme)108 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)106 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme and ABL) 107 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice 108 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 109 109 110 110 LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB … … 113 113 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 114 114 LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB 115 LOGICAL :: ln_tpot !!GS: flag to compute or not potential temperature 115 116 ! 116 117 INTEGER :: nhumi ! choice of the bulk algorithm … … 170 171 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 171 172 & cn_dir , rn_zqt, rn_zu, & 172 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, 173 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, ln_tpot, & 173 174 & ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB 174 175 !!--------------------------------------------------------------------- … … 593 594 !#LB: because AGRIF hates functions that return something else than a scalar, need to 594 595 ! use scalar version of gamma_moist() ... 595 DO_2D_11_11 596 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 597 END_2D 596 IF( ln_tpot ) THEN 597 DO_2D_11_11 598 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 599 END_2D 600 ELSE 601 ztpot = ptair(:,:) 602 ENDIF 598 603 ENDIF 599 604 … … 663 668 ELSE !== BLK formulation ==! turbulent fluxes computation 664 669 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 665 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &666 & wndm(:,:), zU_zu(:,:), pslp(:,:), &667 & taum(:,:), psen(:,:), zqla(:,:), &668 & pEvap=pevp(:,:), prhoa=rhoa(:,:) )670 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 671 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 672 & taum(:,:), psen(:,:), zqla(:,:), & 673 & pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 669 674 670 675 zqla(:,:) = zqla(:,:) * tmask(:,:,1) … … 893 898 894 899 ! local scalars ( place there for vector optimisation purposes) 895 !IF (ln_abl) rhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI)896 900 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 897 901 … … 913 917 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' & 914 918 & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) 915 ELSE 919 ELSE ! ln_abl 916 920 zztmp1 = 11637800.0_wp 917 921 zztmp2 = -5897.8_wp
Note: See TracChangeset
for help on using the changeset viewer.
