Changeset 13540 for NEMO/branches/2020/r12377_ticket2386/src/ABL
- Timestamp:
- 2020-09-29T12:41:06+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/r12377_ticket2386
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/r12377_ticket2386
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 8 9 9 # SETTE 10 ^/utils/CI/sette@ HEADsette10 ^/utils/CI/sette@13507 sette
-
- Property svn:externals
-
NEMO/branches/2020/r12377_ticket2386/src/ABL/abl.F90
r12511 r13540 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/branches/2020/r12377_ticket2386/src/ABL/ablmod.F90
r12511 r13540 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 ! use rn_ ?fac20 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 DO_2D _11_11129 DO_2D( 1, 1, 1, 1 ) 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( 1, 1, 1, 1 ) 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) ) 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) ) 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) ) 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) ) 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 !------------- 431 DO_2D _01_01484 !------------- 485 DO_2D( 0, 1, 0, 1 ) 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 !------------- 455 DO_2D _11_11507 !------------- 508 DO_2D( 1, 1, 1, 1 ) 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. ) 480 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 ! 482 ! first ABL level 526 !------------- 527 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 528 ! ! 6 *** MPI exchanges 529 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 530 ! 531 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1._wp, v_abl(:,:,:,nt_a) , 'T', -1._wp ) 532 CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1._wp , tq_abl(:,:,:,nt_a,jp_qa), 'T', 1._wp , kfillmode = jpfillnothing ) ! ++++ this should not be needed... 533 ! 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 520 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 580 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 581 ! 582 DO_2D( 1, 1, 1, 1 ) 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 531 DO_2D _01_01532 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) )590 591 DO_2D( 0, 1, 0, 1 ) 592 zwnd_i(ji,jj) = u_abl(ji ,jj,2,nt_a) - 0.5_wp * ( pssu(ji ,jj) + pssu(ji-1,jj) ) 593 zwnd_j(ji,jj) = v_abl(ji,jj ,2,nt_a) - 0.5_wp * ( pssv(ji,jj ) + pssv(ji,jj-1) ) 534 594 END_2D 535 ! 536 CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1. , zwnd_j(:,:) , 'T', -1.)595 ! 596 CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1.0_wp, zwnd_j(:,:) , 'T', -1.0_wp ) 537 597 ! 538 598 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 539 DO_2D _11_11599 DO_2D( 1, 1, 1, 1 ) 540 600 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 541 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) )! * msk_abl(ji,jj)601 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 542 602 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 543 603 544 604 pwndm (ji,jj) = zcff 545 605 ptaum (ji,jj) = zztmp * zcff … … 550 610 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 551 611 ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 552 DO_2D _00_00612 DO_2D( 0, 0, 0, 0 ) 553 613 zcff = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) ) 554 614 zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj)) … … 559 619 END_2D 560 620 ! 561 CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1. , ptauj(:,:), 'V', -1.)621 CALL lbc_lnk_multi( 'ablmod', ptaui(:,:), 'U', -1.0_wp, ptauj(:,:), 'V', -1.0_wp ) 562 622 563 623 CALL iom_put( "taum_oce", ptaum ) 564 624 565 625 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: ' )626 CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau : ', mask1=umask, & 627 & tab2d_2=ptauj , clinfo2=' vtau : ', mask2=vmask ) 628 CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm : ' ) 569 629 ENDIF 570 630 571 631 #if defined key_si3 572 ! ------------------------------------------------------------ ! 573 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 574 ! ------------------------------------------------------------ ! 575 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 580 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 581 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 582 & * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 583 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 584 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 585 & * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 586 END_2D 587 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 588 ! 589 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' & 590 & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' ) 632 ! ------------------------------------------------------------ ! 633 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 634 ! ------------------------------------------------------------ ! 635 DO_2D( 0, 0, 0, 0 ) 636 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) & 637 & * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) ) 638 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj) ) & 639 & * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) ) 640 END_2D 641 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp ) 642 ! 643 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: putaui : ' & 644 & , tab2d_2=ptauj_ice , clinfo2=' pvtaui : ' ) 645 ! ------------------------------------------------------------ ! 646 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 647 ! ------------------------------------------------------------ ! 648 DO_2D( 0, 0, 0, 0 ) 649 650 zztmp1 = 0.5_wp * ( u_abl(ji+1,jj ,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 651 zztmp2 = 0.5_wp * ( v_abl(ji ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 652 653 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 654 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 655 & * ( zztmp1 - pssu_ice(ji,jj) ) 656 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 657 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 658 & * ( zztmp2 - pssv_ice(ji,jj) ) 659 END_2D 660 CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice,'V', -1.0_wp ) 661 ! 662 IF(sn_cfctl%l_prtctl) THEN 663 CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask, & 664 & tab2d_2=ptauj_ice , clinfo2=' vtau_ice : ', mask2=vmask ) 665 END IF 591 666 #endif 592 667 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 593 668 ! ! 8 *** Swap time indices for the next timestep 594 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 595 nt_n = 1 + MOD( kt, 2)596 nt_a = 1 + MOD( kt+1, 2)597 ! 669 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 670 nt_n = 1 + MOD( nt_n, 2) 671 nt_a = 1 + MOD( nt_a, 2) 672 ! 598 673 !--------------------------------------------------------------------------------------------------- 599 674 END SUBROUTINE abl_stp 600 675 !=================================================================================================== 601 602 603 604 605 606 607 608 609 610 611 612 613 614 676 615 677 … … 634 696 !! (= Kz dz[Ub] * dz[Un] ) 635 697 !! --------------------------------------------------------------------- 636 INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn 698 INTEGER :: ji, jj, jk, tind, jbak, jkup, jkdwn 637 699 INTEGER, DIMENSION(1:jpi ) :: ikbl 638 700 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, zwndj701 REAL(wp) :: zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2 702 REAL(wp) :: zdU2, zdV2, zbuoy1, zbuoy2 ! zbuoy for BL89 703 REAL(wp) :: zwndi, zwndj 642 704 REAL(wp), DIMENSION(1:jpi, 1:jpka) :: zsh2 643 705 REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: zbn2 644 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF 706 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: zFC, zRH, zCF 645 707 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_a 646 708 REAL(wp), DIMENSION(1:jpi,1:jpka ) :: z_elem_b … … 648 710 LOGICAL :: ln_Patankar = .FALSE. 649 711 LOGICAL :: ln_dumpvar = .FALSE. 650 LOGICAL , DIMENSION(1:jpi ) :: ln_foundl 712 LOGICAL , DIMENSION(1:jpi ) :: ln_foundl 651 713 ! 652 714 tind = nt_n … … 660 722 !------------- 661 723 ! 662 ! Compute vertical shear 724 ! Compute vertical shear 663 725 DO jk = 2, jpkam1 664 DO ji = 1, jpi665 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 726 DO ji = 1, jpi 727 zcff = 1.0_wp / e3w_abl( jk )**2 728 zdU = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 667 729 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 730 zsh2(ji,jk) = zdU+zdV !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 ) 731 END DO 732 END DO 671 733 ! 672 734 ! Compute brunt-vaisala frequency 673 735 DO jk = 2, jpkam1 674 DO ji = 1,jpi 675 zcff = grav * itvref / e3w_abl( jk ) 736 DO ji = 1,jpi 737 zcff = grav * itvref / e3w_abl( jk ) 676 738 zcff1 = tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk , tind, jp_ta) 677 739 zcff2 = tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa) & … … 679 741 zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 ) !<-- zbn2 defined on (2,jpi) 680 742 END DO 681 END DO 743 END DO 682 744 ! 683 745 ! Terms for the tridiagonal problem 684 746 DO jk = 2, jpkam1 685 DO ji = 1, jpi686 zshear = zsh2( ji, jk )! zsh2 is already multiplied by Avm_abl at this point687 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation688 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )689 690 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-diagonal691 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-diagonal747 DO ji = 1, jpi 748 zshear = zsh2( ji, jk ) ! zsh2 is already multiplied by Avm_abl at this point 749 zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk ) ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 750 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 751 752 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 753 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 754 IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE 693 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-side755 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 756 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) ! diagonal 757 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 758 ELSE 697 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 * zbuoy700 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side759 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 760 & + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk) & ! diagonal 761 & - e3w_abl(jk) * rDt_abl * zbuoy 762 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * zshear ) ! right-hand-side 701 763 END IF 702 764 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) 715 zsh2(ji, jpka) = zsh2( ji , jpkam1) 716 END DO 765 END DO 766 767 DO ji = 1,jpi ! vector opt. 768 zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min ) 769 zetop = tke_min 770 771 z_elem_a ( ji, 1 ) = 0._wp 772 z_elem_c ( ji, 1 ) = 0._wp 773 z_elem_b ( ji, 1 ) = 1._wp 774 tke_abl ( ji, jj, 1, nt_a ) = zesrf 775 776 !++ Top Neumann B.C. 777 !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 ) 778 !z_elem_c ( ji, jpka ) = 0._wp 779 !z_elem_b ( ji, jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 780 !tke_abl ( ji, jj, jpka, nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 781 782 !++ Top Dirichlet B.C. 783 z_elem_a ( ji, jpka ) = 0._wp 784 z_elem_c ( ji, jpka ) = 0._wp 785 z_elem_b ( ji, jpka ) = 1._wp 786 tke_abl ( ji, jj, jpka, nt_a ) = zetop 787 788 zbn2 ( ji, jj, 1 ) = zbn2 ( ji, jj, 2 ) 789 zsh2 ( ji, 1 ) = zsh2 ( ji, 2 ) 790 zbn2 ( ji, jj, jpka ) = zbn2 ( ji, jj, jpkam1 ) 791 zsh2 ( ji, jpka ) = zsh2 ( ji , jpkam1 ) 792 END DO 717 793 !! 718 794 !! Matrix inversion … … 720 796 DO ji = 1,jpi 721 797 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 798 zCF (ji, 1 ) = - zcff * z_elem_c( ji, 1 ) 799 tke_abl(ji,jj,1,nt_a) = zcff * tke_abl ( ji, jj, 1, nt_a ) 800 END DO 801 802 DO jk = 2, jpka 727 803 DO ji = 1,jpi 728 804 zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) ) … … 732 808 END DO 733 809 END DO 734 735 DO jk = jpkam1,1,-1 810 811 DO jk = jpkam1,1,-1 736 812 DO ji = 1,jpi 737 813 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 814 END DO 739 815 END DO 740 741 !!FL should not be needed because of Patankar procedure 816 817 !!FL should not be needed because of Patankar procedure 742 818 tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min ) 743 819 … … 745 821 !! Diagnose PBL height 746 822 !! ---------------------------------------------------------- 747 748 749 ! 823 824 825 ! 750 826 ! arrays zRH, zFC and zCF are available at this point 751 827 ! and zFC(:, 1 ) = 0. 752 828 ! diagnose PBL height based on zsh2 and zbn2 753 829 zFC ( : ,1) = 0._wp 754 ikbl( 1:jpi ) = 0 755 830 ikbl( 1:jpi ) = 0 831 756 832 DO jk = 2,jpka 757 DO ji = 1, jpi 833 DO ji = 1, jpi 758 834 zcff = ghw_abl( jk-1 ) 759 835 zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) ) … … 781 857 ELSE 782 858 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 ... 859 END IF 860 END DO 861 !------------- 862 END DO 863 !------------- 864 ! 865 ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 790 866 IF(ln_smth_pblh) THEN 791 CALL lbc_lnk( 'ablmod', pblh, 'T', 1. )867 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing) 792 868 CALL smooth_pblh( pblh, msk_abl ) 793 CALL lbc_lnk( 'ablmod', pblh, 'T', 1. )869 CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing) 794 870 ENDIF 795 871 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 799 875 SELECT CASE ( nn_amxl ) 800 876 ! 801 CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom 802 # define zlup zRH 803 # define zldw zFC 877 CASE ( 0 ) ! Deardroff 80 length-scale bounded by the distance to surface and bottom 878 # define zlup zRH 879 # define zldw zFC 804 880 DO jj = 1, jpj ! outer loop 805 881 ! 806 882 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 883 mxld_abl( ji, jj, 1 ) = mxl_min 884 mxld_abl( ji, jj, jpka ) = mxl_min 885 mxlm_abl( ji, jj, 1 ) = mxl_min 886 mxlm_abl( ji, jj, jpka ) = mxl_min 887 zldw ( ji, 1 ) = zrough(ji,jj) * rn_Lsfc 888 zlup ( ji, jpka ) = mxl_min 889 END DO 890 ! 891 DO jk = 2, jpkam1 892 DO ji = 1, jpi 893 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 894 mxlm_abl( ji, jj, jk ) = MAX( mxl_min, & 895 & SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) ) 896 END DO 897 END DO 820 898 ! 821 899 ! 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 900 DO jk = jpkam1,1,-1 901 DO ji = 1, jpi 902 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 903 END DO 904 END DO 827 905 ! 828 906 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 907 DO ji = 1, jpi 908 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 909 END DO 910 END DO 911 ! 912 ! DO jk = 1, jpka 913 ! DO ji = 1, jpi 914 ! mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 915 ! mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ), zlup( ji, jk ) ) 916 ! END DO 917 ! END DO 833 918 ! 834 919 DO jk = 1, jpka 835 920 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 921 ! zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 922 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 923 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 924 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 925 END DO 926 END DO 927 ! 928 END DO 929 # undef zlup 930 # undef zldw 931 ! 932 ! 933 CASE ( 1 ) ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom 934 # define zlup zRH 935 # define zldw zFC 936 DO jj = 1, jpj ! outer loop 937 ! 938 DO jk = 2, jpkam1 939 DO ji = 1,jpi 940 zcff = 1.0_wp / e3w_abl( jk )**2 941 zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 942 zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 943 zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 944 ENDDO 945 ENDDO 946 ! 947 DO ji = 1, jpi 948 zcff = zrough(ji,jj) * rn_Lsfc 949 mxld_abl ( ji, jj, 1 ) = zcff 950 mxld_abl ( ji, jj, jpka ) = mxl_min 951 mxlm_abl ( ji, jj, 1 ) = zcff 952 mxlm_abl ( ji, jj, jpka ) = mxl_min 953 zldw ( ji, 1 ) = zcff 954 zlup ( ji, jpka ) = mxl_min 955 END DO 956 ! 957 DO jk = 2, jpkam1 958 DO ji = 1, jpi 959 zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 960 zcff = 2.0_wp*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) & 961 & + SQRT(rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.0_wp*zbuoy ) ) 962 mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff ) 963 END DO 964 END DO 965 ! 966 ! Limit mxl 967 DO jk = jpkam1,1,-1 968 DO ji = 1, jpi 969 zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 970 END DO 971 END DO 972 ! 973 DO jk = 2, jpka 974 DO ji = 1, jpi 975 zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 976 END DO 977 END DO 978 ! 979 DO jk = 1, jpka 980 DO ji = 1, jpi 981 !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 982 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 983 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 984 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 985 !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) ) 986 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 987 END DO 988 END DO 989 ! 990 END DO 991 # undef zlup 992 # undef zldw 863 993 ! 864 994 CASE ( 2 ) ! Bougeault & Lacarrere 89 length-scale 865 995 ! 866 # define zlup zRH 867 # define zldw zFC 996 # define zlup zRH 997 # define zldw zFC 868 998 ! zCF is used for matrix inversion 869 ! 999 ! 870 1000 DO jj = 1, jpj ! outer loop 871 872 DO ji = 1, jpi 873 zlup( ji, 1 ) = mxl_min 874 zldw( ji, 1 ) = mxl_min 1001 1002 DO ji = 1, jpi 1003 zcff = zrough(ji,jj) * rn_Lsfc 1004 zlup( ji, 1 ) = zcff 1005 zldw( ji, 1 ) = zcff 875 1006 zlup( ji, jpka ) = mxl_min 876 zldw( ji, jpka ) = mxl_min 877 END DO 878 1007 zldw( ji, jpka ) = mxl_min 1008 END DO 1009 879 1010 DO jk = 2,jpka-1 880 1011 DO ji = 1, jpi 881 1012 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 1013 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 1014 END DO 1015 END DO 885 1016 !! 886 1017 !! BL89 search for lup 887 !! ---------------------------------------------------------- 888 DO jk=2,jpka-1 1018 !! ---------------------------------------------------------- 1019 DO jk=2,jpka-1 889 1020 ! 890 1021 DO ji = 1, jpi … … 892 1023 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 893 1024 ln_foundl(ji ) = .false. 894 END DO 895 ! 1025 END DO 1026 ! 896 1027 DO jkup=jk+1,jpka-1 897 1028 DO ji = 1, jpi 1029 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1030 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 898 1031 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)) )1032 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & 1033 & + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) 901 1034 IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 902 1035 zcff2 = ghw_abl(jkup ) - ghw_abl(jk) 903 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 1036 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 904 1037 zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / & 905 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 906 zlup(ji,jk) = zcff 1038 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 1039 zlup(ji,jk) = zcff 1040 zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk) 907 1041 ln_foundl(ji) = .true. 908 1042 END IF … … 910 1044 END DO 911 1045 ! 912 END DO 1046 END DO 913 1047 !! 914 1048 !! BL89 search for ldwn 915 !! ---------------------------------------------------------- 916 DO jk=2,jpka-1 1049 !! ---------------------------------------------------------- 1050 DO jk=2,jpka-1 917 1051 ! 918 1052 DO ji = 1, jpi … … 920 1054 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 921 1055 ln_foundl(ji ) = .false. 922 END DO 923 ! 1056 END DO 1057 ! 924 1058 DO jkdwn=jk-1,1,-1 925 DO ji = 1, jpi 1059 DO ji = 1, jpi 1060 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1061 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 926 1062 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 1063 & * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 1064 + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) 1065 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 930 1066 zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 931 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 1067 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 932 1068 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 ! 1069 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1070 zldw(ji,jk) = zcff 1071 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1072 ln_foundl(ji) = .true. 1073 END IF 1074 END DO 1075 END DO 1076 ! 940 1077 END DO 941 1078 942 1079 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 1080 DO ji = 1, jpi 1081 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1082 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1083 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1084 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1085 END DO 1086 END DO 947 1087 948 1088 END DO 949 # undef zlup 950 # undef zldw 951 ! 952 END SELECT 1089 # undef zlup 1090 # undef zldw 1091 ! 1092 CASE ( 3 ) ! Bougeault & Lacarrere 89 length-scale 1093 ! 1094 # define zlup zRH 1095 # define zldw zFC 1096 ! zCF is used for matrix inversion 1097 ! 1098 DO jj = 1, jpj ! outer loop 1099 ! 1100 DO jk = 2, jpkam1 1101 DO ji = 1,jpi 1102 zcff = 1.0_wp / e3w_abl( jk )**2 1103 zdU = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk , tind) )**2 1104 zdV = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk , tind) )**2 1105 zsh2(ji,jk) = SQRT(zdU+zdV) !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 1106 ENDDO 1107 ENDDO 1108 zsh2(:, 1) = zsh2( :, 2) 1109 zsh2(:, jpka) = zsh2( :, jpkam1) 1110 1111 DO ji = 1, jpi 1112 zcff = zrough(ji,jj) * rn_Lsfc 1113 zlup( ji, 1 ) = zcff 1114 zldw( ji, 1 ) = zcff 1115 zlup( ji, jpka ) = mxl_min 1116 zldw( ji, jpka ) = mxl_min 1117 END DO 1118 1119 DO jk = 2,jpka-1 1120 DO ji = 1, jpi 1121 zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 1122 zldw(ji,jk) = ghw_abl(jk ) - ghw_abl( 1) 1123 END DO 1124 END DO 1125 !! 1126 !! BL89 search for lup 1127 !! ---------------------------------------------------------- 1128 DO jk=2,jpka-1 1129 ! 1130 DO ji = 1, jpi 1131 zCF(ji,1:jpka) = 0._wp 1132 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 1133 ln_foundl(ji ) = .false. 1134 END DO 1135 ! 1136 DO jkup=jk+1,jpka-1 1137 DO ji = 1, jpi 1138 zbuoy1 = MAX( zbn2(ji,jj,jkup ), rsmall ) 1139 zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 1140 zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 1141 & ( zbuoy1*(ghw_abl(jkup )-ghw_abl(jk)) & 1142 & + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) & 1143 & + 0.5_wp * e3t_abl(jkup) * rn_Rod * & 1144 & ( SQRT(tke_abl( ji, jj, jkup , nt_a ))*zsh2(ji,jkup ) & 1145 & + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) ) 1146 1147 IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1148 zcff2 = ghw_abl(jkup ) - ghw_abl(jk) 1149 zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 1150 zcff = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) / & 1151 & ( zCF(ji,jkup) - zCF(ji,jkup-1) ) 1152 zlup(ji,jk) = zcff 1153 zlup(ji,jk) = ghw_abl(jkup ) - ghw_abl(jk) 1154 ln_foundl(ji) = .true. 1155 END IF 1156 END DO 1157 END DO 1158 ! 1159 END DO 1160 !! 1161 !! BL89 search for ldwn 1162 !! ---------------------------------------------------------- 1163 DO jk=2,jpka-1 1164 ! 1165 DO ji = 1, jpi 1166 zCF(ji,1:jpka) = 0._wp 1167 zCF(ji, jk ) = - tke_abl( ji, jj, jk, nt_a ) 1168 ln_foundl(ji ) = .false. 1169 END DO 1170 ! 1171 DO jkdwn=jk-1,1,-1 1172 DO ji = 1, jpi 1173 zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 1174 zbuoy2 = MAX( zbn2(ji,jj,jkdwn ), rsmall ) 1175 zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1) & 1176 & * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 1177 & +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn )) ) & 1178 & + 0.5_wp * e3t_abl(jkup) * rn_Rod * & 1179 & ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) & 1180 & + SQRT(tke_abl( ji, jj, jkdwn , nt_a ))*zsh2(ji,jkdwn ) ) 1181 1182 IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 1183 zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 1184 zcff1 = ghw_abl(jk) - ghw_abl(jkdwn ) 1185 zcff = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) / & 1186 & ( zCF(ji,jkdwn+1) - zCF(ji,jkdwn) ) 1187 zldw(ji,jk) = zcff 1188 zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn ) 1189 ln_foundl(ji) = .true. 1190 END IF 1191 END DO 1192 END DO 1193 ! 1194 END DO 1195 1196 DO jk = 1, jpka 1197 DO ji = 1, jpi 1198 !zcff = 2.*SQRT(2.)*( zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp) )**(-3._wp/2._wp) 1199 zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 1200 mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 1201 mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ), zlup( ji, jk ) ), mxl_min ) 1202 END DO 1203 END DO 1204 1205 END DO 1206 # undef zlup 1207 # undef zldw 1208 ! 1209 ! 1210 END SELECT 953 1211 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 954 1212 ! ! Finalize the computation of turbulent visc./diff. 955 1213 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 956 1214 957 1215 !------------- 958 1216 DO jj = 1, jpj ! outer loop 959 1217 !------------- 960 DO jk = 1, jpka 1218 DO jk = 1, jpka 961 1219 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 1220 zcff = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk ) & 1221 & * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) ) 1222 zcff2 = 1. / ( 1. + zcff ) !<-- phi_z(z) 1223 zcff = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) ) 1224 !!FL: MAX function probably useless because of the definition of mxl_min 967 1225 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 1226 Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak ) 1227 END DO 1228 END DO 1229 !------------- 1230 END DO 973 1231 !------------- 974 1232 … … 988 1246 !! 989 1247 !! --------------------------------------------------------------------- 990 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk 991 1248 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: msk 1249 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d 992 1250 INTEGER :: ji,jj 993 994 995 1251 REAL(wp) :: smth_a, smth_b 1252 REAL(wp), DIMENSION(jpi,jpj) :: zdX,zdY,zFX,zFY 1253 REAL(wp) :: zumsk,zvmsk 996 1254 !! 997 1255 !!========================================================= … … 1001 1259 smth_b = 1._wp / 4._wp 1002 1260 ! 1003 DO_2D _11_101261 DO_2D( 1, 1, 1, 0 ) 1004 1262 zumsk = msk(ji,jj) * msk(ji+1,jj) 1005 1263 zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji ,jj ) ) * zumsk 1006 1264 END_2D 1007 1008 DO_2D_10_11 1265 1266 DO_2D( 1, 0, 1, 1 ) 1009 1267 zvmsk = msk(ji,jj) * msk(ji,jj+1) 1010 1268 zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji ,jj ) ) * zvmsk 1011 1012 1013 DO_2D_10_00 1269 END_2D 1270 1271 DO_2D( 1, 0, 0, 0 ) 1014 1272 zFY ( ji, jj ) = zdY ( ji, jj ) & 1015 1273 & + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) & 1016 1274 & - (zdX ( ji, jj ) - zdX( ji-1, jj )) ) 1017 1018 1019 DO_2D _00_101275 END_2D 1276 1277 DO_2D( 0, 0, 1, 0 ) 1020 1278 zFX( ji, jj ) = zdX( ji, jj ) & 1021 1279 & + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) & … … 1023 1281 END_2D 1024 1282 1025 DO_2D _00_001283 DO_2D( 0, 0, 0, 0 ) 1026 1284 pvar2d( ji ,jj ) = pvar2d( ji ,jj ) & 1027 1285 & + msk(ji,jj) * smth_b * ( & … … 1029 1287 & +zFY( ji, jj ) - zFY( ji, jj-1 ) ) 1030 1288 END_2D 1031 !! 1289 1032 1290 !--------------------------------------------------------------------------------------------------- 1033 1291 END SUBROUTINE smooth_pblh -
NEMO/branches/2020/r12377_ticket2386/src/ABL/ablrst.F90
r11945 r13540 74 74 ENDIF 75 75 ! 76 CALL iom_open( TRIM(clpath)//TRIM(clname), numraw, ldwrt = .TRUE., kdlev = jpka )76 CALL iom_open( TRIM(clpath)//TRIM(clname), numraw, ldwrt = .TRUE., kdlev = jpka, cdcomp = 'ABL' ) 77 77 lrst_abl = .TRUE. 78 78 ENDIF … … 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 ! … … 146 146 ENDIF 147 147 148 CALL iom_open ( TRIM(cn_ablrst_indir)//'/'//cn_ablrst_in, numrar , kdlev = jpka)148 CALL iom_open ( TRIM(cn_ablrst_indir)//'/'//cn_ablrst_in, numrar ) 149 149 150 150 ! Time info … … 165 165 166 166 ! --- mandatory fields --- ! 167 CALL iom_get( numrar, jpdom_auto glo, 'u_abl', u_abl(:,:,:,nt_n ))168 CALL iom_get( numrar, jpdom_auto glo, 'v_abl', v_abl(:,:,:,nt_n ))169 CALL iom_get( numrar, jpdom_auto glo, 't_abl', tq_abl(:,:,:,nt_n,jp_ta) )170 CALL iom_get( numrar, jpdom_auto glo, 'q_abl', tq_abl(:,:,:,nt_n,jp_qa) )171 CALL iom_get( numrar, jpdom_auto glo, 'tke_abl', tke_abl(:,:,:,nt_n ) )172 CALL iom_get( numrar, jpdom_auto glo, 'avm_abl', avm_abl(:,:,: ) )173 CALL iom_get( numrar, jpdom_auto glo, 'avt_abl', avt_abl(:,:,: ) )174 CALL iom_get( numrar, jpdom_auto glo, 'mxl_abl', mxl_abl(:,:,: ) )175 CALL iom_get( numrar, jpdom_auto glo, 'pblh', pblh(:,: ) )167 CALL iom_get( numrar, jpdom_auto, 'u_abl', u_abl(:,:,:,nt_n ), cd_type = 'U', psgn = -1._wp ) 168 CALL iom_get( numrar, jpdom_auto, 'v_abl', v_abl(:,:,:,nt_n ), cd_type = 'V', psgn = -1._wp ) 169 CALL iom_get( numrar, jpdom_auto, 't_abl', tq_abl(:,:,:,nt_n,jp_ta) ) 170 CALL iom_get( numrar, jpdom_auto, 'q_abl', tq_abl(:,:,:,nt_n,jp_qa) ) 171 CALL iom_get( numrar, jpdom_auto, 'tke_abl', tke_abl(:,:,:,nt_n ) ) 172 CALL iom_get( numrar, jpdom_auto, 'avm_abl', avm_abl(:,:,: ) ) 173 CALL iom_get( numrar, jpdom_auto, 'avt_abl', avt_abl(:,:,: ) ) 174 CALL iom_get( numrar, jpdom_auto,'mxld_abl',mxld_abl(:,:,: ) ) 175 CALL iom_get( numrar, jpdom_auto, 'pblh', pblh(:,: ) ) 176 176 CALL iom_delay_rst( 'READ', 'ABL', numrar ) ! read only abl delayed global communication variables 177 177 -
NEMO/branches/2020/r12377_ticket2386/src/ABL/par_abl.F90
r12511 r13540 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 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/r12377_ticket2386/src/ABL/sbcabl.F90
r12511 r13540 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 REWIND( numnam_ref )! 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 80 ! 81 REWIND( numnam_cfg )! Namelist namsbc_abl in configuration namelist : ABL parameters81 ! Namelist namsbc_abl in configuration namelist : ABL parameters 82 82 READ ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 ) 83 83 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist' ) … … 166 166 rn_Sch = rn_ce / rn_cm 167 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 168 170 169 171 IF(lwp) THEN … … 172 174 WRITE(numout,*) ' ~~~~~~~~~~~' 173 175 IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 174 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 ' 175 179 WRITE(numout,*) ' Minimum value of atmospheric TKE = ',tke_min,' m^2 s^-2' 176 180 WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m' … … 179 183 WRITE(numout,*) ' Constant for Schmidt number = ',rn_Sch 180 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 181 187 END IF 182 188 … … 203 209 ! ABL timestep 204 210 rDt_abl = nn_fsbc * rn_Dt 211 IF(lwp) WRITE(numout,*) ' ABL timestep = ', rDt_abl,' s' 205 212 206 213 ! Check parameters for dynamics … … 249 256 zcff = 2._wp * omega * SIN( rad * 90._wp ) !++ fmax 250 257 rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 251 !!GS: alternative shape252 !rest_eq(:,:) = SIN( 0.5_wp*rpi*(zcff - ABS(ff_t(:,:))) / (zcff - 3.e-5) )**8253 !WHERE(ABS(ff_t(:,:)).LE.3.e-5) rest_eq(:,:) = 1._wp254 258 ELSE 255 259 rest_eq(:,:) = 1._wp … … 264 268 265 269 ! Initialize the time index for now time (nt_n) and after time (nt_a) 266 nt_n = 1 + MOD( nit000 , 2) 267 nt_a = 1 + MOD( nit000+1, 2) 270 nt_n = 1; nt_a = 2 268 271 269 272 ! initialize ABL from data or restart 270 IF( ln_rstart ) THEN273 IF( ln_rstart_abl ) THEN 271 274 CALL abl_rst_read 272 275 ELSE 273 276 CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 274 277 275 u_abl(:,:,:,nt_n ) = sf(jp_wndi)%fnow(:,:,:)276 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(:,:,:) 277 280 tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 278 281 tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) … … 281 284 avm_abl(:,:,: ) = avm_bak 282 285 avt_abl(:,:,: ) = avt_bak 283 mxl_abl(:,:,: ) = mxl_min284 286 pblh (:,: ) = ghw_abl( 3 ) !<-- assume that the pbl contains 3 grid points 285 287 u_abl (:,:,:,nt_a ) = 0._wp … … 287 289 tq_abl (:,:,:,nt_a,: ) = 0._wp 288 290 tke_abl(:,:,:,nt_a ) = 0._wp 291 292 mxlm_abl(:,:,: ) = mxl_min 293 mxld_abl(:,:,: ) = mxl_min 289 294 ENDIF 290 291 rhoa(:,:) = rho_air( tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), sf(jp_slp)%fnow(:,:,1) ) !!GS: rhoa must be (re)computed here here to avoid division by zero in blk_ice_1 (TBI)292 295 293 296 END SUBROUTINE sbc_abl_init … … 330 333 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 331 334 332 !!------------------------------------------------------------------------------------------- 333 !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields 334 !!------------------------------------------------------------------------------------------- 335 336 CALL blk_oce_1( kt, u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in 337 & tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in 338 & sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in 339 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in 340 & tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out 341 342 #if defined key_si3 343 CALL blk_ice_1( u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in 344 & tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in 345 & sf(jp_slp)%fnow(:,:,1) , u_ice, v_ice, tm_su , & ! <<= in 346 & pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui ) ! <<= out 347 #endif 348 349 !!------------------------------------------------------------------------------------------- 350 !! 3 - Advance ABL variables from now (n) to after (n+1) 351 !!------------------------------------------------------------------------------------------- 352 353 CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq, & ! <<= in 354 & sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:), & ! <<= in 355 & sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:), & ! <<= in 356 & sf(jp_slp )%fnow(:,:,1), & ! <<= in 357 & sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:), & ! <<= in 358 & zcd_du, zsen, zevp, & ! <=> in/out 359 & wndm, utau, vtau, taum & ! =>> out 360 #if defined key_si3 361 & , tm_su, u_ice, v_ice, zssqi, zcd_dui & ! <<= in 362 & , zseni, zevpi, wndm_ice, ato_i & ! <<= in 363 & , utau_ice, vtau_ice & ! =>> out 364 #endif 365 & ) 366 !!------------------------------------------------------------------------------------------- 367 !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since 368 !! time swap is done in abl_stp 369 !!------------------------------------------------------------------------------------------- 370 371 CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta), & 372 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1), & 373 & sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1), & 374 & tsk_m, zsen, zevp ) 375 376 CALL abl_rst_opn( kt ) ! Open abl restart file (if necessary) 377 IF( lrst_abl ) CALL abl_rst_write( kt ) ! -- abl restart file 378 379 #if defined key_si3 380 ! Avoid a USE abl in icesbc module 381 sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta); sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa) 382 #endif 335 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 336 337 !!------------------------------------------------------------------------------------------- 338 !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields 339 !!------------------------------------------------------------------------------------------- 340 341 CALL blk_oce_1( kt, u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in 342 & tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in 343 & sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in 344 & sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in 345 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in 346 & tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out 347 348 #if defined key_si3 349 CALL blk_ice_1( u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in 350 & tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in 351 & sf(jp_slp)%fnow(:,:,1) , u_ice, v_ice, tm_su , & ! <<= in 352 & pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui ) ! <<= out 353 #endif 354 355 !!------------------------------------------------------------------------------------------- 356 !! 3 - Advance ABL variables from now (n) to after (n+1) 357 !!------------------------------------------------------------------------------------------- 358 359 CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq, & ! <<= in 360 & sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:), & ! <<= in 361 & sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:), & ! <<= in 362 & sf(jp_slp )%fnow(:,:,1), & ! <<= in 363 & sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:), & ! <<= in 364 & zcd_du, zsen, zevp, & ! <=> in/out 365 & wndm, utau, vtau, taum & ! =>> out 366 #if defined key_si3 367 & , tm_su, u_ice, v_ice, zssqi, zcd_dui & ! <<= in 368 & , zseni, zevpi, wndm_ice, ato_i & ! <<= in 369 & , utau_ice, vtau_ice & ! =>> out 370 #endif 371 & ) 372 !!------------------------------------------------------------------------------------------- 373 !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since 374 !! time swap is done in abl_stp 375 !!------------------------------------------------------------------------------------------- 376 377 CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta), & 378 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1), & 379 & sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1), & 380 & tsk_m, zsen, zevp ) 381 382 CALL abl_rst_opn( kt ) ! Open abl restart file (if necessary) 383 IF( lrst_abl ) CALL abl_rst_write( kt ) ! -- abl restart file 384 385 #if defined key_si3 386 ! Avoid a USE abl in icesbc module 387 sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta); sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa) 388 #endif 389 END IF 383 390 384 391 END SUBROUTINE sbc_abl
Note: See TracChangeset
for help on using the changeset viewer.