[13826] | 1 | MODULE zdfmfc |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE zdfmfc *** |
---|
| 4 | !! Ocean physics: Mass-Flux scheme parameterization of Convection: |
---|
| 5 | !! Non-local transport for the convective ocean boundary |
---|
| 6 | !! layer. Subgrid-scale large eddies are represented by a |
---|
| 7 | !! mass-flux contribution (ln_zdfmfc = .TRUE.) |
---|
| 8 | !!====================================================================== |
---|
| 9 | !! History : NEMO ! |
---|
| 10 | !! 3.6 ! 2016-06 (H. Giordani, R. Bourdallé-Badie) Original code |
---|
| 11 | !! 4.2 ! 2020-12 (H. Giordani, R. Bourdallé-Badie) adapt to NEM04.2 |
---|
| 12 | !!---------------------------------------------------------------------- |
---|
| 13 | !!---------------------------------------------------------------------- |
---|
| 14 | !! tra_mfc : Compute the Mass Flux and trends of T/S |
---|
| 15 | !! diag_mfc : Modify diagonal of trazdf Matrix |
---|
| 16 | !! rhs_mfc : Modify RHS of trazdf Matrix |
---|
| 17 | !! zdf_mfc_init : initialization, namelist read, and parameters control |
---|
| 18 | !!---------------------------------------------------------------------- |
---|
| 19 | ! |
---|
| 20 | USE oce ! ocean dynamics and active tracers |
---|
| 21 | USE dom_oce ! ocean space and time domain |
---|
| 22 | USE domvvl ! ocean space and time domain : variable volume layer |
---|
| 23 | USE domzgr |
---|
| 24 | USE zdf_oce ! ocean vertical physics |
---|
| 25 | USE sbc_oce ! surface boundary condition: ocean |
---|
| 26 | USE phycst ! physical constants |
---|
| 27 | USE eosbn2 ! equation of state (eos routine) |
---|
| 28 | USE zdfmxl ! mixed layer |
---|
| 29 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
| 30 | USE lib_mpp ! MPP manager |
---|
| 31 | USE prtctl ! Print control |
---|
| 32 | USE in_out_manager ! I/O manager |
---|
| 33 | USE iom ! I/O manager library |
---|
| 34 | USE timing ! Timing |
---|
| 35 | USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) |
---|
| 36 | |
---|
| 37 | IMPLICIT NONE |
---|
| 38 | PRIVATE |
---|
| 39 | |
---|
| 40 | PUBLIC tra_mfc ! routine called in step module |
---|
| 41 | PUBLIC diag_mfc ! routine called in trazdf module |
---|
| 42 | PUBLIC rhs_mfc ! routine called in trazdf module |
---|
| 43 | PUBLIC zdf_mfc_init ! routine called in nemo module |
---|
| 44 | ! |
---|
| 45 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: edmfa, edmfb, edmfc !: diagonal term of the matrix. |
---|
| 46 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: edmftra !: y term for matrix inversion |
---|
| 47 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: edmfm !: y term for matrix inversion |
---|
| 48 | ! |
---|
| 49 | !! ** Namelist namzdf_edmf ** |
---|
| 50 | REAL(wp) :: rn_cemf ! entrain of T/S |
---|
| 51 | REAL(wp) :: rn_cwmf ! detrain of T/S |
---|
| 52 | REAL(wp) :: rn_cent ! entrain of the convective mass flux |
---|
| 53 | REAL(wp) :: rn_cdet ! detrain of the convective mass flux |
---|
| 54 | REAL(wp) :: rn_cap ! Factor of computation for convective area (negative => area constant) |
---|
| 55 | REAL(wp) :: App_max ! Maximum of the convective area |
---|
| 56 | LOGICAL, PUBLIC, SAVE :: ln_edmfuv !: EDMF flag for velocity ! |
---|
| 57 | ! |
---|
| 58 | !! * Substitutions |
---|
| 59 | # include "do_loop_substitute.h90" |
---|
| 60 | # include "domzgr_substitute.h90" |
---|
| 61 | !!---------------------------------------------------------------------- |
---|
| 62 | !! NEMO/OCE 4.2 , NEMO Consortium (2018) |
---|
| 63 | !! $Id: zdfmfc.F90 13783 2020-20-02 15:30:22Z rbourdal $ |
---|
| 64 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
| 65 | !!---------------------------------------------------------------------- |
---|
| 66 | CONTAINS |
---|
| 67 | |
---|
| 68 | INTEGER FUNCTION zdf_mfc_alloc() |
---|
| 69 | !!---------------------------------------------------------------------- |
---|
| 70 | !! *** FUNCTION zdf_edmf_alloc *** |
---|
| 71 | !!---------------------------------------------------------------------- |
---|
| 72 | ALLOCATE( edmfa(jpi,jpj,jpk) , edmfb(jpi,jpj,jpk) , edmfc(jpi,jpj,jpk) & |
---|
| 73 | & , edmftra(jpi,jpj,jpk,2), edmfm(jpi,jpj,jpk) , STAT= zdf_mfc_alloc ) |
---|
| 74 | ! |
---|
| 75 | IF( lk_mpp ) CALL mpp_sum ( 'zdfmfc', zdf_mfc_alloc ) |
---|
| 76 | IF( zdf_mfc_alloc /= 0 ) CALL ctl_warn('zdf_mfc_alloc: failed to allocate arrays') |
---|
| 77 | END FUNCTION zdf_mfc_alloc |
---|
| 78 | |
---|
| 79 | |
---|
| 80 | SUBROUTINE tra_mfc( kt, Kmm, pts, Krhs ) |
---|
| 81 | !!---------------------------------------------------------------------- |
---|
| 82 | !! *** ROUTINE zdf_mfc *** |
---|
| 83 | !! |
---|
| 84 | !! ** Purpose : Compute a mass flux, depending on surface flux, over |
---|
| 85 | !! the instable part of the water column. |
---|
| 86 | !! |
---|
| 87 | !! ** Method : Compute surface instability and mix tracers until stable level |
---|
| 88 | !! |
---|
| 89 | !! |
---|
| 90 | !! ** Action : Compute convection plume and (ta,sa)-trends for trazdf (EDMF scheme) |
---|
| 91 | !! |
---|
| 92 | !! References : |
---|
| 93 | !! Giordani, Bourdallé-Badie and Madec JAMES 2020 |
---|
| 94 | !!---------------------------------------------------------------------- |
---|
| 95 | !!---------------------------------------------------------------------- |
---|
| 96 | INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices |
---|
| 97 | REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation |
---|
| 98 | REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: ztsp ! T/S of the plume |
---|
| 99 | REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: ztse ! T/S at W point |
---|
| 100 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp ! |
---|
| 101 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp2 ! |
---|
| 102 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zapp ! |
---|
| 103 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zedmf ! |
---|
| 104 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zepsT, zepsW ! |
---|
| 105 | ! |
---|
| 106 | REAL(wp), DIMENSION(jpi,jpj) :: zustar, zustar2 ! |
---|
| 107 | REAL(wp), DIMENSION(jpi,jpj) :: zuws, zvws, zsws, zfnet ! |
---|
| 108 | REAL(wp), DIMENSION(jpi,jpj) :: zfbuo, zrautbm1, zrautb, zraupl |
---|
| 109 | REAL(wp), DIMENSION(jpi,jpj) :: zwpsurf ! |
---|
| 110 | REAL(wp), DIMENSION(jpi,jpj) :: zop0 , zsp0 ! |
---|
| 111 | REAL(wp), DIMENSION(jpi,jpj) :: zrwp_0, zrwp2_0 ! |
---|
| 112 | REAL(wp), DIMENSION(jpi,jpj) :: zapp0 ! |
---|
| 113 | REAL(wp), DIMENSION(jpi,jpj) :: zphp, zph, zphpm1, zphm1, zNHydro |
---|
| 114 | REAL(wp), DIMENSION(jpi,jpj) :: zhcmo ! |
---|
| 115 | ! |
---|
| 116 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zn2 ! N^2 |
---|
| 117 | REAL(wp), DIMENSION(jpi,jpj,2 ) :: zab, zabm1, zabp ! alpha and beta |
---|
| 118 | |
---|
[13832] | 119 | REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value |
---|
| 120 | |
---|
[13826] | 121 | REAL(wp) :: zrho, zrhop |
---|
| 122 | REAL(wp) :: zcnh, znum, zden, zcoef1, zcoef2 |
---|
| 123 | REAL(wp) :: zca, zcb, zcd, zrw, zxl, zcdet, zctre |
---|
| 124 | REAL(wp) :: zaw, zbw, zxw |
---|
| 125 | REAL(wp) :: alpha |
---|
| 126 | ! |
---|
| 127 | INTEGER, INTENT(in ) :: kt ! ocean time-step index ! |
---|
| 128 | ! |
---|
| 129 | INTEGER :: ji, jj, jk ! dummy loop arguments |
---|
| 130 | ! |
---|
| 131 | !------------------------------------------------------------------ |
---|
| 132 | ! Initialisation of coefficients |
---|
| 133 | !------------------------------------------------------------------ |
---|
| 134 | zca = 1._wp |
---|
| 135 | zcb = 1._wp |
---|
| 136 | zcd = 1._wp |
---|
| 137 | |
---|
| 138 | !------------------------------------------------------------------ |
---|
| 139 | ! Surface boundary condition |
---|
| 140 | !------------------------------------------------------------------ |
---|
| 141 | ! surface Stress |
---|
| 142 | !-------------------- |
---|
| 143 | zuws(:,:) = utau(:,:) * r1_rho0 |
---|
| 144 | zvws(:,:) = vtau(:,:) * r1_rho0 |
---|
| 145 | zustar2(:,:) = SQRT(zuws(:,:)*zuws(:,:)+zvws(:,:)*zvws(:,:)) |
---|
| 146 | zustar(:,:) = SQRT(zustar2(:,:)) |
---|
| 147 | |
---|
| 148 | ! Heat Flux |
---|
| 149 | !-------------------- |
---|
| 150 | zfnet(:,:) = qns(:,:) + qsr(:,:) |
---|
| 151 | zfnet(:,:) = zfnet(:,:) / (rho0 * rcp) |
---|
| 152 | |
---|
| 153 | ! Water Flux |
---|
| 154 | !--------------------- |
---|
| 155 | zsws(:,:) = emp(:,:) |
---|
| 156 | |
---|
| 157 | !------------------------------------------- |
---|
| 158 | ! Initialisation of prognostic variables |
---|
| 159 | !------------------------------------------- |
---|
[13888] | 160 | zrwp (:,:,:) = 0._wp ; zrwp2(:,:,:) = 0._wp ; zedmf(:,:,:) = 0._wp |
---|
| 161 | zph (:,:) = 0._wp ; zphm1(:,:) = 0._wp ; zphpm1(:,:) = 0._wp |
---|
| 162 | ztsp(:,:,:,:)= 0._wp |
---|
[13826] | 163 | |
---|
| 164 | ! Tracers inside plume (ztsp) and environment (ztse) |
---|
| 165 | ztsp(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1) |
---|
| 166 | ztsp(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) |
---|
| 167 | ztse(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1) |
---|
| 168 | ztse(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) |
---|
| 169 | |
---|
| 170 | CALL eos( ztse(:,:,1,:) , zrautb(:,:) ) |
---|
| 171 | CALL eos( ztsp(:,:,1,:) , zraupl(:,:) ) |
---|
| 172 | |
---|
| 173 | !------------------------------------------- |
---|
| 174 | ! Boundary Condition of Mass Flux (plume velo.; convective area, entrain/detrain) |
---|
| 175 | !------------------------------------------- |
---|
| 176 | zhcmo(:,:) = e3t(:,:,1,Kmm) |
---|
[13888] | 177 | zfbuo(:,:) = 0._wp |
---|
[13826] | 178 | WHERE ( ABS(zrautb(:,:)) > 1.e-20 ) zfbuo(:,:) = & |
---|
| 179 | & grav * ( 2.e-4_wp *zfnet(:,:) - 7.6E-4_wp*pts(:,:,1,jp_sal,Kmm)*zsws(:,:)/zrautb(:,:)) * zhcmo(:,:) |
---|
| 180 | |
---|
| 181 | zedmf(:,:,1) = -0.065_wp*(ABS(zfbuo(:,:)))**(1._wp/3._wp)*SIGN(1.,zfbuo(:,:)) |
---|
| 182 | zedmf(:,:,1) = MAX(0., zedmf(:,:,1)) |
---|
| 183 | |
---|
| 184 | zwpsurf(:,:) = 2._wp/3._wp*zustar(:,:) + 2._wp/3._wp*ABS(zfbuo(:,:))**(1._wp/3._wp) |
---|
| 185 | zwpsurf(:,:) = MAX(1.e-5_wp,zwpsurf(:,:)) |
---|
| 186 | zwpsurf(:,:) = MIN(1.,zwpsurf(:,:)) |
---|
| 187 | |
---|
| 188 | zapp(:,:,:) = App_max |
---|
| 189 | WHERE(zwpsurf .NE. 0.) zapp(:,:,1) = MIN(MAX(0.,zedmf(:,:,1)/zwpsurf(:,:)), App_max) |
---|
| 190 | |
---|
[13955] | 191 | zedmf(:,:,1) = 0._wp |
---|
| 192 | zrwp (:,:,1) = 0._wp |
---|
| 193 | zrwp2(:,:,1) = 0._wp |
---|
[13888] | 194 | zepsT(:,:,:) = 0.001_wp |
---|
| 195 | zepsW(:,:,:) = 0.001_wp |
---|
[13826] | 196 | |
---|
| 197 | |
---|
| 198 | !-------------------------------------------------------------- |
---|
| 199 | ! Compute plume properties |
---|
| 200 | ! In the same loop on vert. levels computation of: |
---|
| 201 | ! - Vertical velocity: zWp |
---|
| 202 | ! - Convective Area: zAp |
---|
| 203 | ! - Tracers properties inside the plume (if necessary): ztp |
---|
| 204 | !--------------------------------------------------------------- |
---|
| 205 | |
---|
| 206 | DO jk= 2, jpk |
---|
| 207 | |
---|
| 208 | ! Compute the buoyancy acceleration on T-points at jk-1 |
---|
| 209 | zrautbm1(:,:) = zrautb(:,:) |
---|
| 210 | CALL eos( pts (:,:,jk ,:,Kmm) , zrautb(:,:) ) |
---|
| 211 | CALL eos( ztsp(:,:,jk-1,: ) , zraupl(:,:) ) |
---|
| 212 | |
---|
| 213 | zphm1(:,:) = zphm1(:,:) + grav * zrautbm1(:,:) * e3t(:,:,jk-1, Kmm) |
---|
| 214 | zphpm1(:,:) = zphpm1(:,:) + grav * zraupl(:,:) * e3t(:,:,jk-1, Kmm) |
---|
| 215 | zph(:,:) = zphm1(:,:) + grav * zrautb(:,:) * e3t(:,:,jk , Kmm) |
---|
[13832] | 216 | zph(:,:) = MAX( zph(:,:), zepsilon) |
---|
[13826] | 217 | |
---|
| 218 | WHERE(zrautbm1 .NE. 0.) zfbuo(:,:) = grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:) |
---|
| 219 | |
---|
| 220 | DO_2D( 0, 0, 0, 0 ) |
---|
| 221 | |
---|
| 222 | ! Compute Environment of Plume. Interpolation T/S (before time step) on W-points |
---|
| 223 | zrw = (gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm)) & |
---|
| 224 | & / (gdept(ji,jj,jk,Kmm) - gdept(ji,jj,jk-1,Kmm)) |
---|
| 225 | ztse(ji,jj,jk,:) = (pts(ji,jj,jk,:,Kmm) * zrw + pts(ji,jj,jk-1,:,Kmm)*(1._wp - zrw) )*tmask(ji,jj,jk) |
---|
| 226 | |
---|
| 227 | !--------------------------------------------------------------- |
---|
| 228 | ! Compute the vertical velocity on W-points |
---|
| 229 | !--------------------------------------------------------------- |
---|
| 230 | |
---|
| 231 | ! Non-hydrostatic pressure terms in the wp2 equation |
---|
[13888] | 232 | zcnh = 0.2_wp |
---|
[13826] | 233 | znum = 0.5_wp + zcnh - & |
---|
| 234 | (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) & |
---|
| 235 | *e3t(ji,jj,jk-1,Kmm)*0.5_wp |
---|
| 236 | zden = 0.5_wp + zcnh + & |
---|
| 237 | (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) & |
---|
| 238 | *e3t(ji,jj,jk-1,Kmm)*0.5_wp |
---|
| 239 | |
---|
| 240 | zcoef1 = zca*e3t(ji,jj,jk-1,Kmm) / zden |
---|
| 241 | zcoef2 = znum/zden |
---|
| 242 | |
---|
| 243 | ! compute wp2 |
---|
| 244 | zrwp2(ji,jj,jk) = zcoef1*zfbuo(ji,jj) & |
---|
| 245 | + zcoef2*zrwp2(ji,jj,jk-1) |
---|
| 246 | zrwp2(ji,jj,jk) = MAX ( zrwp2(ji,jj,jk)*wmask(ji,jj,jk) , 0.) |
---|
| 247 | zrwp (ji,jj,jk) = SQRT( zrwp2(ji,jj,jk) ) |
---|
| 248 | |
---|
| 249 | !---------------------------------------------------------------------------------- |
---|
| 250 | ! Compute convective area on W-point |
---|
| 251 | ! Compute vertical profil of the convective area with mass conservation hypothesis |
---|
| 252 | ! If rn_cap negative => constant value on the water column. |
---|
| 253 | !---------------------------------------------------------------------------------- |
---|
| 254 | IF( rn_cap .GT. 0. ) THEN |
---|
| 255 | |
---|
| 256 | zxw = MAX(zrwp(ji,jj,jk-1), zrwp(ji,jj,jk) ) |
---|
| 257 | IF( zxw > 0. ) THEN |
---|
| 258 | |
---|
| 259 | zxl = (zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk))/(e3t(ji,jj,jk-1,Kmm)*zxw) |
---|
| 260 | IF (zxl .LT. 0._wp) THEN |
---|
| 261 | zctre = -1.*rn_cap*zxl |
---|
[13888] | 262 | zcdet = 0._wp |
---|
[13826] | 263 | ELSE |
---|
[13888] | 264 | zctre = 0._wp |
---|
[13826] | 265 | zcdet = rn_cap*zxl |
---|
| 266 | END IF |
---|
| 267 | zapp(ji,jj,jk) = zapp(ji,jj,jk-1)* & |
---|
| 268 | & (1._wp + (zxl + zctre - zcdet )*e3t(ji,jj,jk-1,Kmm)) |
---|
| 269 | ELSE |
---|
| 270 | zapp(ji,jj,jk) = App_max |
---|
| 271 | END IF |
---|
| 272 | zapp(ji,jj,jk) = MIN( MAX(zapp(ji,jj,jk),0.), App_max) |
---|
| 273 | ELSE |
---|
| 274 | zapp(ji,jj,jk) = -1. * rn_cap |
---|
| 275 | END IF |
---|
| 276 | |
---|
| 277 | ! Compute Mass Flux on W-point |
---|
| 278 | zedmf(ji,jj,jk) = -zapp(ji,jj,jk) * zrwp(ji,jj,jk)* wmask(ji,jj,jk) |
---|
| 279 | |
---|
| 280 | ! Compute Entrainment coefficient |
---|
| 281 | IF(rn_cemf .GT. 0.) THEN |
---|
[13888] | 282 | zxw = 0.5_wp*(zrwp(ji,jj,jk-1)+ zrwp(ji,jj,jk) ) |
---|
| 283 | zepsT(ji,jj,jk) = 0.01_wp |
---|
[13826] | 284 | IF( zxw > 0. ) THEN |
---|
| 285 | zepsT(ji,jj,jk) = zepsT(ji,jj,jk) + & |
---|
[13888] | 286 | & ABS( zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk) ) & |
---|
[13826] | 287 | & / ( e3t(ji,jj,jk-1,Kmm) * zxw ) |
---|
| 288 | zepsT(ji,jj,jk) = zepsT(ji,jj,jk) * rn_cemf * wmask(ji,jj,jk) |
---|
| 289 | ENDIF |
---|
| 290 | ELSE |
---|
| 291 | zepsT(ji,jj,jk) = -rn_cemf |
---|
| 292 | ENDIF |
---|
| 293 | |
---|
| 294 | ! Compute the detrend coef for velocity (on W-point and not T-points, bug ???) |
---|
| 295 | IF(rn_cwmf .GT. 0.) THEN |
---|
| 296 | zepsW(ji,jj,jk) = rn_cwmf * zepsT(ji,jj,jk) |
---|
| 297 | ELSE |
---|
| 298 | zepsW(ji,jj,jk) = -rn_cwmf |
---|
| 299 | ENDIF |
---|
| 300 | |
---|
| 301 | !--------------------------------------------------------------- |
---|
| 302 | ! Compute the plume properties on T-points |
---|
| 303 | !--------------------------------------------------------------- |
---|
[13888] | 304 | IF(zrwp (ji,jj,jk) .LT. 1.e-12_wp .AND. zrwp (ji,jj,jk-1) .LT. 1.e-12_wp) THEN |
---|
[13826] | 305 | ztsp(ji,jj,jk-1,jp_tem) = pts(ji,jj,jk-1,jp_tem,Kmm) |
---|
| 306 | ztsp(ji,jj,jk-1,jp_sal) = pts(ji,jj,jk-1,jp_sal,Kmm) |
---|
| 307 | ENDIF |
---|
| 308 | |
---|
| 309 | zcoef1 = (1._wp-zepsT(ji,jj,jk)*(1._wp-zrw)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk ) ) & |
---|
| 310 | & / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk) ) |
---|
| 311 | ! |
---|
| 312 | zcoef2 = zepsT(ji,jj,jk)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk) & |
---|
| 313 | & / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk)) |
---|
| 314 | ! |
---|
| 315 | ztsp(ji,jj,jk,jp_tem) = (zcoef1 * ztsp(ji,jj,jk-1,jp_tem) + & |
---|
| 316 | & zcoef2 * ztse(ji,jj,jk ,jp_tem) )*tmask(ji,jj,jk) |
---|
| 317 | ztsp(ji,jj,jk,jp_sal) = (zcoef1 * ztsp(ji,jj,jk-1,jp_sal) + & |
---|
| 318 | & zcoef2 * ztse(ji,jj,jk ,jp_sal) )*tmask(ji,jj,jk) |
---|
| 319 | |
---|
| 320 | END_2D |
---|
| 321 | END DO ! end of loop on jpk |
---|
| 322 | |
---|
| 323 | ! Compute Mass Flux on T-point |
---|
| 324 | DO jk=1,jpk-1 |
---|
[13888] | 325 | edmfm(:,:,jk) = (zedmf(:,:,jk+1) + zedmf(:,:,jk) )*0.5_wp |
---|
[13826] | 326 | END DO |
---|
[13832] | 327 | edmfm(:,:,jpk) = zedmf(:,:,jpk) |
---|
[13826] | 328 | |
---|
| 329 | ! Save variable (on T point) |
---|
| 330 | CALL iom_put( "mf_Tp" , ztsp(:,:,:,jp_tem) ) ! Save plume temperature |
---|
| 331 | CALL iom_put( "mf_Sp" , ztsp(:,:,:,jp_sal) ) ! Save plume salinity |
---|
| 332 | CALL iom_put( "mf_mf" , edmfm(:,:,:) ) ! Save Mass Flux |
---|
| 333 | ! Save variable (on W point) |
---|
| 334 | CALL iom_put( "mf_wp" , zrwp (:,:,:) ) ! Save convective velocity in the plume |
---|
| 335 | CALL iom_put( "mf_app", zapp (:,:,:) ) ! Save convective area |
---|
| 336 | |
---|
| 337 | !================================================================================= |
---|
| 338 | ! Computation of a tridiagonal matrix and right hand side terms of the linear system |
---|
| 339 | !================================================================================= |
---|
[13888] | 340 | edmfa(:,:,:) = 0._wp |
---|
| 341 | edmfb(:,:,:) = 0._wp |
---|
| 342 | edmfc(:,:,:) = 0._wp |
---|
| 343 | edmftra(:,:,:,:) = 0._wp |
---|
[13826] | 344 | |
---|
| 345 | !--------------------------------------------------------------- |
---|
| 346 | ! Diagonal terms |
---|
| 347 | !--------------------------------------------------------------- |
---|
| 348 | DO jk=1,jpk-1 |
---|
[13888] | 349 | edmfa(:,:,jk) = 0._wp |
---|
[13826] | 350 | edmfb(:,:,jk) = -edmfm(:,:,jk ) / e3w(:,:,jk+1,Kmm) |
---|
| 351 | edmfc(:,:,jk) = edmfm(:,:,jk+1) / e3w(:,:,jk+1,Kmm) |
---|
| 352 | END DO |
---|
[13832] | 353 | edmfa(:,:,jpk) = -edmfm(:,:,jpk-1) / e3w(:,:,jpk,Kmm) |
---|
| 354 | edmfb(:,:,jpk) = edmfm(:,:,jpk ) / e3w(:,:,jpk,Kmm) |
---|
[13888] | 355 | edmfc(:,:,jpk) = 0._wp |
---|
[13826] | 356 | |
---|
| 357 | !--------------------------------------------------------------- |
---|
| 358 | ! right hand side term for Temperature |
---|
| 359 | !--------------------------------------------------------------- |
---|
| 360 | DO jk=1,jpk-1 |
---|
| 361 | edmftra(:,:,jk,1) = - edmfm(:,:,jk ) * ztsp(:,:,jk ,jp_tem) / e3w(:,:,jk+1,Kmm) & |
---|
| 362 | & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_tem) / e3w(:,:,jk+1,Kmm) |
---|
| 363 | END DO |
---|
| 364 | edmftra(:,:,jpk,1) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_tem) / e3w(:,:,jpk,Kmm) & |
---|
| 365 | & + edmfm(:,:,jpk ) * ztsp(:,:,jpk ,jp_tem) / e3w(:,:,jpk,Kmm) |
---|
| 366 | |
---|
| 367 | !--------------------------------------------------------------- |
---|
| 368 | ! Right hand side term for Salinity |
---|
| 369 | !--------------------------------------------------------------- |
---|
| 370 | DO jk=1,jpk-1 |
---|
| 371 | edmftra(:,:,jk,2) = - edmfm(:,:,jk ) * ztsp(:,:,jk ,jp_sal) / e3w(:,:,jk+1,Kmm) & |
---|
| 372 | & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_sal) / e3w(:,:,jk+1,Kmm) |
---|
| 373 | END DO |
---|
| 374 | edmftra(:,:,jpk,2) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_sal) / e3w(:,:,jpk,Kmm) & |
---|
| 375 | & + edmfm(:,:,jpk ) * ztsp(:,:,jpk ,jp_sal) / e3w(:,:,jpk,Kmm) |
---|
| 376 | ! |
---|
| 377 | ! |
---|
| 378 | CALL lbc_lnk_multi( 'zdfmfc', edmfm,'T',1., edmfa,'T',1., edmfb,'T',1., edmfc,'T',1., edmftra(:,:,:,jp_tem),'T',1., edmftra(:,:,:,jp_sal),'T',1.) |
---|
| 379 | ! |
---|
| 380 | END SUBROUTINE tra_mfc |
---|
| 381 | |
---|
| 382 | |
---|
| 383 | SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa ) |
---|
| 384 | |
---|
| 385 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: zdiagi, zdiagd, zdiags ! inout: tridaig. terms |
---|
| 386 | REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step |
---|
| 387 | INTEGER , INTENT(in ) :: Kaa ! ocean time level indices |
---|
| 388 | |
---|
| 389 | INTEGER :: ji, jj, jk ! dummy loop arguments |
---|
| 390 | |
---|
| 391 | DO_3D( 0, 0, 0, 0, 1, jpkm1 ) |
---|
| 392 | zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk) |
---|
| 393 | zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk) |
---|
| 394 | zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk) |
---|
| 395 | END_3D |
---|
| 396 | |
---|
| 397 | END SUBROUTINE diag_mfc |
---|
| 398 | |
---|
| 399 | SUBROUTINE rhs_mfc( zrhs, jjn ) |
---|
| 400 | |
---|
| 401 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: zrhs ! inout: rhs trend |
---|
| 402 | INTEGER , INTENT(in ) :: jjn ! tracer indices |
---|
| 403 | |
---|
| 404 | INTEGER :: ji, jj, jk ! dummy loop arguments |
---|
| 405 | |
---|
| 406 | DO_3D( 0, 0, 0, 0, 1, jpkm1 ) |
---|
| 407 | zrhs(ji,jj,jk) = zrhs(ji,jj,jk) + edmftra(ji,jj,jk,jjn) |
---|
| 408 | END_3D |
---|
| 409 | |
---|
| 410 | END SUBROUTINE rhs_mfc |
---|
| 411 | |
---|
| 412 | |
---|
| 413 | |
---|
| 414 | SUBROUTINE zdf_mfc_init |
---|
| 415 | !!---------------------------------------------------------------------- |
---|
| 416 | !! *** ROUTINE zdf_mfc_init *** |
---|
| 417 | !! |
---|
| 418 | !! ** Purpose : Initialization of the vertical eddy diffivity and |
---|
| 419 | !! mass flux |
---|
| 420 | !! |
---|
| 421 | !! ** Method : Read the namzdf_mfc namelist and check the parameters |
---|
| 422 | !! called at the first timestep (nit000) |
---|
| 423 | !! |
---|
| 424 | !! ** input : Namlist namzdf_mfc |
---|
| 425 | !! |
---|
| 426 | !! ** Action : Increase by 1 the nstop flag is setting problem encounter |
---|
| 427 | !! |
---|
| 428 | !!---------------------------------------------------------------------- |
---|
| 429 | ! |
---|
| 430 | INTEGER :: jk ! dummy loop indices |
---|
| 431 | INTEGER :: ios ! Local integer output status for namelist read |
---|
| 432 | REAL(wp):: zcr ! local scalar |
---|
| 433 | !! |
---|
| 434 | NAMELIST/namzdf_mfc/ ln_edmfuv, rn_cemf, rn_cwmf, rn_cent, rn_cdet, rn_cap, App_max |
---|
| 435 | !!---------------------------------------------------------- |
---|
| 436 | ! |
---|
| 437 | ! |
---|
| 438 | ! REWIND( numnam_ref ) ! Namelist namzdf_mfc in reference namelist : Vertical eddy diffivity mass flux |
---|
| 439 | READ ( numnam_ref, namzdf_mfc, IOSTAT = ios, ERR = 901) |
---|
| 440 | 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in reference namelist' ) |
---|
| 441 | |
---|
| 442 | ! REWIND( numnam_cfg ) ! Namelist namzdf_mfc in configuration namelist : Vertical eddy diffivity mass flux |
---|
| 443 | READ ( numnam_cfg, namzdf_mfc, IOSTAT = ios, ERR = 902 ) |
---|
| 444 | 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in configuration namelist' ) |
---|
| 445 | IF(lwm) WRITE ( numond, namzdf_mfc ) |
---|
| 446 | |
---|
| 447 | IF(lwp) THEN !* Control print |
---|
| 448 | WRITE(numout,*) |
---|
| 449 | WRITE(numout,*) 'zdf_mfc_init' |
---|
| 450 | WRITE(numout,*) '~~~~~~~~~~~~~' |
---|
| 451 | WRITE(numout,*) ' Namelist namzdf_mfc : set eddy diffusivity Mass Flux Convection' |
---|
| 452 | WRITE(numout,*) ' Apply mass flux on velocities (Not yet avail.) ln_edmfuv = ', ln_edmfuv |
---|
| 453 | WRITE(numout,*) ' Coeff for entrain/detrain T/S of plume (Neg => cte) rn_cemf = ', rn_cemf |
---|
| 454 | WRITE(numout,*) ' Coeff for entrain/detrain Wp of plume (Neg => cte) rn_cwmf = ', rn_cwmf |
---|
| 455 | WRITE(numout,*) ' Coeff for entrain/detrain area of plume rn_cap = ', rn_cap |
---|
| 456 | WRITE(numout,*) ' Coeff for entrain area of plume rn_cent = ', rn_cent |
---|
| 457 | WRITE(numout,*) ' Coeff for detrain area of plume rn_cdet = ', rn_cdet |
---|
| 458 | WRITE(numout,*) ' Max convective area App_max = ', App_max |
---|
| 459 | ENDIF |
---|
| 460 | !* allocate edmf arrays |
---|
| 461 | IF( zdf_mfc_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_edmf_init : unable to allocate arrays' ) |
---|
[13888] | 462 | edmfa(:,:,:) = 0._wp |
---|
| 463 | edmfb(:,:,:) = 0._wp |
---|
| 464 | edmfc(:,:,:) = 0._wp |
---|
| 465 | edmftra(:,:,:,:) = 0._wp |
---|
[13826] | 466 | ! |
---|
| 467 | END SUBROUTINE zdf_mfc_init |
---|
| 468 | |
---|
| 469 | !!====================================================================== |
---|
| 470 | |
---|
| 471 | !!====================================================================== |
---|
| 472 | END MODULE zdfmfc |
---|
| 473 | |
---|
| 474 | |
---|
| 475 | |
---|