MITgcm |
|
|
|
|||
File indexing completed on 2018-03-02 18:37:46 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTCb2ea1d2979 Jean*0001 MODULE SHALLOW_CONV_MOD 0002 0003 !======================================================================= 0004 ! --- SHALLOW CONVECTION MODULE - GFDL SPECTRAL MODEL VERSION 0005 !======================================================================= 0006 0007 use simple_Sat_Vapor_Pres_Mod, ONLY: ESCOMP, DESCOMP 0008 ! use fms_Mod, ONLY: FILE_EXIST, ERROR_MESG, FATAL, & 0009 ! CHECK_NML_ERROR, OPEN_namelist_FILE, & 0010 ! CLOSE_FILE, mpp_pe, mpp_root_pe, stdlog, & 0011 ! write_version_number 0012 0013 use constants_mod, only: Hlv, Cp_air, RDgas, RVgas, Kappa, grav 0014 0015 !--------------------------------------------------------------------- 0016 implicit none 0017 private 0018 !--------------------------------------------------------------------- 0019 0020 public :: SHALLOW_CONV, SHALLOW_CONV_INIT 0021 public :: MYLCL 0022 0023 !--------------------------------------------------------------------- 0024 0025 character(len=128) :: version = '$Id: shallow_conv_mod.F90,v 1.1 2013/05/08 22:14:14 jmc Exp $' 0026 character(len=128) :: tag = '$Name: $' 0027 0028 logical :: do_init = .true. 0029 0030 !--------------------------------------------------------------------- 0031 ! --- CONSTANTS 0032 !--------------------------------------------------------------------- 0033 0034 real :: Hlv_by_Cp, Cp_by_RDgas, omkappa, dovcns, d622, d378 0035 real :: crtkons 0036 real, parameter :: p00 = 1000.0E2 0037 integer :: kctopm1, kctopm2 0038 0039 real, allocatable, dimension(:) :: rhcrit, rhmax, rhmin, delrhc 0040 0041 !--------------------------------------------------------------------- 0042 ! --- NAMELIST 0043 !--------------------------------------------------------------------- 0044 0045 logical :: lipps = .false. 0046 logical :: ldetran = .true. 0047 real :: theqvcr = 0.0 0048 real :: pshalow = 750.0E2 0049 real :: akhsc0 = 5.0 0050 real :: crthum = 0.85 0051 real :: hc = 1.0 0052 integer :: kctop = 3 0053 0054 0055 NAMELIST / shallow_conv_nml / & 0056 lipps, ldetran, theqvcr, pshalow, akhsc0, kctop, crthum, hc 0057 0058 !--------------------------------------------------------------------- 0059 0060 contains 0061 0062 !####################################################################### 0063 !####################################################################### 0064 0065 SUBROUTINE SHALLOW_CONV_INIT( kx ) 0066 0067 !======================================================================= 0068 ! ***** INITIALIZE SHALLOW CONVECTION 0069 !======================================================================= 0070 !--------------------------------------------------------------------- 0071 ! Arguments (Intent in) 0072 ! kx - Number of levels in vertical 0073 !--------------------------------------------------------------------- 0074 integer, intent(in) :: kx 0075 0076 !--------------------------------------------------------------------- 0077 ! (Intent local) 0078 !--------------------------------------------------------------------- 0079 integer :: unit, io, ierr 0080 0081 !===================================================================== 0082 0083 if (.not.do_init) return 0084 0085 !--------------------------------------------------------------------- 0086 ! --- Read namelist 0087 !--------------------------------------------------------------------- 0088 0089 !del if( FILE_EXIST( 'input.nml' ) ) then 0090 ! ------------------------------------- 0091 !del unit = OPEN_namelist_FILE () 0092 !del ierr = 1 0093 !del do while( ierr .ne. 0 ) 0094 !del READ ( unit, nml = shallow_conv_nml, iostat = io, end = 10 ) 0095 !del ierr = CHECK_NML_ERROR(io,'shallow_conv_nml') 0096 !del end do 0097 10 continue 0098 !del CALL CLOSE_FILE ( unit ) 0099 ! ------------------------------------- 0100 !del end if 0101 0102 !--------------------------------------------------------------------- 0103 ! --- Output version 0104 !--------------------------------------------------------------------- 0105 0106 !del call write_version_number(version, tag) 0107 !del if ( mpp_pe() == mpp_root_pe() ) write (stdlog(), nml=shallow_conv_nml) 0108 0109 !--------------------------------------------------------------------- 0110 ! --- Initialize constants 0111 !--------------------------------------------------------------------- 0112 0113 d622 = RDgas / RVgas 0114 d378 = 1.0 - d622 0115 Hlv_by_Cp = Hlv / Cp_air 0116 Cp_by_RDgas = Cp_air / RDgas 0117 omkappa = 1.0 - Kappa 0118 dovcns = -0.5 * grav / RDgas 0119 0120 !---------------------------------------------------- 0121 0122 crtkons = -1.0 * theqvcr * RDgas / grav 0123 if( lipps ) crtkons = 0.0 0124 0125 !---------------------------------------------------- 0126 0127 kctopm1 = kctop - 1 0128 kctopm2 = kctop - 2 0129 0130 !---------------------------------------------------- 0131 0132 allocate( rhcrit(kx) ) 0133 allocate( rhmax(kx) ) 0134 allocate( rhmin(kx) ) 0135 allocate( delrhc(kx) ) 0136 0137 rhcrit = crthum 0138 rhmax = hc 0139 rhmin = 2.0 * rhcrit - hc - 1.0e-15 0140 delrhc = 1.0 / ( rhmax - rhmin ) 0141 0142 !------------------------------------------------------------------- 0143 do_init = .false. 0144 !--------------------------------------------------------------------- 0145 0146 !===================================================================== 0147 end SUBROUTINE SHALLOW_CONV_INIT 0148 0149 !####################################################################### 0150 0151 SUBROUTINE SHALLOW_CONV( Temp, qmix0, pfull, phalf, akhsc, kbot ) 0152 0153 !======================================================================= 0154 ! --- SHALLOW CONVECTION 0155 !======================================================================= 0156 !---------------------------------------------------------------------- 0157 ! Arguments (Intent in) 0158 ! Temp - Temperature 0159 ! qmix0 - Specific humidity 0160 ! pfull - Pressure at full levels 0161 ! phalf - Pressure at half levels 0162 ! kbot - OPTIONAL; lowest model level index (integer) 0163 !---------------------------------------------------------------------- 0164 real, intent(in), dimension(:,:,:) :: Temp, qmix0, pfull, phalf 0165 0166 integer, intent(in), OPTIONAL, dimension(:,:) :: kbot 0167 0168 !---------------------------------------------------------------------- 0169 ! Arguments (Intent out) 0170 ! akhsc - mixing coefficient for heat and moisture 0171 ! due to shallow convection 0172 !---------------------------------------------------------------------- 0173 real, intent(out), dimension(:,:,:) :: akhsc 0174 0175 !---------------------------------------------------------------------- 0176 ! --- local 0177 !---------------------------------------------------------------------- 0178 real, dimension(SIZE(Temp,1),SIZE(Temp,2)) :: & 0179 plcl, rhumjmp, rhumscl, xy1, xy2, xy3 0180 0181 integer, dimension(SIZE(Temp,1),SIZE(Temp,2)) :: & 0182 ksiglcl 0183 0184 real, dimension(SIZE(Temp,1),SIZE(Temp,2),SIZE(Temp,3)) :: & 0185 qmix, dphalf, qsat, rhum, theta, thetav, buoy, xyz1 0186 0187 integer, dimension(SIZE(Temp,1),SIZE(Temp,2),SIZE(Temp,3)) :: & 0188 kbuoy 0189 0190 integer :: k, kx, kxm, kxp, i, ix, j, jx 0191 0192 !======================================================================= 0193 !======================================================================= 0194 0195 ix = SIZE(Temp,1) 0196 jx = SIZE(Temp,2) 0197 kx = SIZE(Temp,3) 0198 0199 kxm = kx - 1 0200 kxp = kx + 1 0201 0202 !======================================================================= 0203 ! --- MOISTURE VARIABLES 0204 !======================================================================= 0205 0206 qmix = qmix0 0207 qmix = MAX( qmix, 1.0E-6 ) 0208 qmix = MIN( qmix, 0.2 ) 0209 0210 ! --- saturation mixing ratio 0211 CALL ESCOMP( Temp, qsat ) 0212 xyz1 = pfull - d378 * qsat 0213 qsat = d622 * qsat / xyz1 0214 0215 ! --- relative humidity 0216 rhum = qmix / qsat 0217 0218 !======================================================================= 0219 ! --- POTENTIAL TEMPERATURE 0220 !======================================================================= 0221 0222 theta = Temp * ( ( p00 / pfull )**Kappa ) 0223 0224 !======================================================================= 0225 ! --- CALCULATE THE LIFTING CONDENSATION LEVEL, IE CLOUB BASE 0226 !======================================================================= 0227 0228 if( PRESENT( kbot ) ) then 0229 do j = 1,jx 0230 do i = 1,ix 0231 k = kbot(i,j) 0232 xy1(i,j) = Temp(i,j,k) 0233 xy2(i,j) = MIN( qmix(i,j,k), qsat(i,j,k) ) 0234 xy3(i,j) = pfull(i,j,k) 0235 end do 0236 end do 0237 else 0238 xy1(:,:) = Temp(:,:,kx) 0239 xy2(:,:) = MIN( qmix(:,:,kx), qsat(:,:,kx) ) 0240 xy3(:,:) = pfull(:,:,kx) 0241 end if 0242 0243 CALL MYLCL( xy1, xy2, xy3, phalf, plcl, ksiglcl ) 0244 0245 !======================================================================= 0246 ! --- INITALIZE 0247 !======================================================================= 0248 0249 kbuoy = kxp 0250 akhsc = 0.0 0251 0252 !======================================================================= 0253 ! --- BUOYANCY 0254 !======================================================================= 0255 0256 !--------------------------------------------------------------------- 0257 ! --- DEFAULT: 0258 ! --- BASED ON EQUIVALENT POTENTIAL TEMPERATURE GRADIENT 0259 !--------------------------------------------------------------------- 0260 0261 if( .not. lipps ) then 0262 ! %%%%%%%%%%%%%%%%%%%%%%% 0263 0264 ! --- Vertical differential of pressure 0265 dphalf(:,:,1:kxm) = pfull(:,:,2:kx) - pfull(:,:,1:kxm) 0266 0267 if( PRESENT( kbot ) ) then 0268 dphalf(:,:,1:kxm) = MAX( dphalf(:,:,1:kxm), 1.0e-5 ) 0269 end if 0270 0271 ! --- Equivalent potential temperature 0272 xyz1 = ( Hlv_by_Cp * qmix ) / Temp 0273 ! thetav = theta * (1.0 + xyz1 ) 0274 thetav = theta * EXP( xyz1 ) 0275 0276 ! --- Equivalent potential temperature gradient 0277 xyz1(:,:,2:kx) = thetav(:,:,2:kx) - thetav(:,:,1:kxm) 0278 xyz1(:,:,2:kx) = xyz1(:,:,2:kx) / dphalf(:,:,1:kxm) 0279 buoy(:,:,2:kxm) = 0.5*( xyz1(:,:,2:kxm) + xyz1(:,:,3:kx) ) 0280 0281 ! %%%%%%%%%%%%%%%%%%%%%%% 0282 endif 0283 0284 !--------------------------------------------------------------------- 0285 ! --- OPTION: 0286 ! --- BUOYANCY ALA FRANK LIPPS 0287 !--------------------------------------------------------------------- 0288 0289 if( lipps ) then 0290 ! %%%%%%%%%%%%%%%%%%%%%%% 0291 0292 ! --- Virtual potential temperature gardient 0293 thetav = theta * ( 1.0 + 0.608 * qmix ) 0294 0295 ! --- Buoyancy 0296 do k = kctopm2,kxm 0297 ! ------------------- 0298 xy1(:,:) = Hlv / ( Cp_air*RVgas*Temp(:,:,k)*Temp(:,:,k) ) + 1.0 / qsat(:,:,k) 0299 0300 xy2(:,:) = ( grav / ( Cp_air*Temp(:,:,k) ) ) * & 0301 ( Hlv / ( RVgas*Temp(:,:,k) ) - Cp_by_RDgas ) / xy1(:,:) 0302 0303 xy3(:,:) = ( thetav(:,:,k+1) - thetav(:,:,k ) ) / phalf(:,:,k+1) + & 0304 ( thetav(:,:,k ) - thetav(:,:,k-1) ) / phalf(:,:,k ) 0305 0306 buoy(:,:,k) = ( hlv / ( cp_air*Temp(:,:,k) ) - 1.608 ) * xy2(:,:) 0307 buoy(:,:,k) = buoy(:,:,k) - dovcns*pfull(:,:,k) * xy3(:,:) / Temp(:,:,k) 0308 ! ------------------- 0309 end do 0310 0311 ! %%%%%%%%%%%%%%%%%%%%%%% 0312 endif 0313 0314 !======================================================================= 0315 ! --- COMPUTE THE LEVEL OF NO BUOYANCY 0316 ! --- RETAIN ONLY THE LOWEST CONTIGUOUS BUOYANT SHALLOW CONVECTIVE LAYER. 0317 !======================================================================= 0318 0319 do k = kctopm1,kxm 0320 ! ------------------- 0321 where ( ( pfull(:,:,k) >= pshalow ) .and. & 0322 ( pfull(:,:,k) <= plcl(:,:) ) .and. & 0323 ( buoy(:,:,k) >= crtkons ) ) 0324 kbuoy(:,:,k) = k 0325 endwhere 0326 ! ------------------- 0327 end do 0328 0329 do k = kctopm1,kxm 0330 ! ------------------- 0331 where( ( pfull(:,:,k) < plcl(:,:) ) .and. & 0332 ( kbuoy(:,:,k) == kxp ) .and. & 0333 ( kbuoy(:,:,k-1) == k-1 ) ) 0334 kbuoy(:,:,k-1) = kxp 0335 endwhere 0336 ! ------------------- 0337 end do 0338 0339 !======================================================================= 0340 ! --- SHALLOW CONVECTION WILL OCCUR AT LEVELS WHERE KBUOY <= KSIGLCL 0341 !======================================================================= 0342 0343 do k = kctopm1,kxm 0344 ! ------------------- 0345 where( kbuoy(:,:,k) <= ksiglcl(:,:) ) 0346 akhsc(:,:,k+1) = akhsc0 0347 endwhere 0348 ! ------------------- 0349 end do 0350 0351 !======================================================================= 0352 ! --- DETRAINMENT THRU INVERSION LAYER 0353 !======================================================================= 0354 0355 !--------------------------------------------------------------------- 0356 ! --- DEFAULT: 0357 ! --- ENHANCED DETRAINMENT THRU INVERSION LAYER. 0358 !--------------------------------------------------------------------- 0359 0360 if( ldetran) then 0361 ! %%%%%%%%%%%%%%%%%%%%%%%% 0362 0363 do k = kctopm1,kxm 0364 ! ------------------- 0365 where( ( kbuoy(:,:,k) == k ) .and. & 0366 ( kbuoy(:,:,k-1) == kxp ) .and. & 0367 ( pfull(:,:,k) >= pshalow ) ) 0368 akhsc(:,:,k) = 0.2 * akhsc0 0369 akhsc(:,:,k+1) = 0.6 * akhsc0 0370 endwhere 0371 ! ------------------- 0372 end do 0373 0374 do k = kctopm1,kxm 0375 ! ------------------- 0376 where( ( pfull(:,:,k) <= plcl(:,:) ) .and. & 0377 ( pfull(:,:,k+1) > plcl(:,:) ) .and. & 0378 ( kbuoy(:,:,k) == k ) ) 0379 akhsc(:,:,k+1) = 0.2 * akhsc0 0380 endwhere 0381 ! ------------------- 0382 end do 0383 0384 ! %%%%%%%%%%%%%%%%%%%%%%%% 0385 endif 0386 0387 !--------------------------------------------------------------------- 0388 ! --- OPTION: 0389 ! --- NORMAL DETRAINMENT THRU INVERSION LAYER 0390 !--------------------------------------------------------------------- 0391 0392 if( .not. ldetran ) then 0393 ! %%%%%%%%%%%%%%%%%%%%%%%% 0394 0395 rhumscl = 0.0 0396 rhumjmp = 0.0 0397 0398 do k = kctopm1,kxm 0399 ! ------------------- 0400 where( ( kbuoy(:,:,k) == k ) .and. & 0401 ( buoy(:,:,k-1) < crtkons ) .and. & 0402 ( buoy(:,:,k) >= crtkons ) ) 0403 rhumjmp(:,:) = rhum(:,:,k) - rhum(:,:,k-1) 0404 rhumscl(:,:) = delrhc(k) * ( rhum(:,:,k) - rhmin(k) ) 0405 endwhere 0406 ! ------------------- 0407 end do 0408 0409 rhumscl = MIN( 1.0, rhumscl ) 0410 rhumscl = MAX( 0.0, rhumscl ) 0411 rhumjmp = MIN( 1.0, rhumjmp ) 0412 rhumjmp = MAX( 0.0, rhumjmp ) 0413 0414 do k = kctopm1,kxm 0415 ! ------------------- 0416 where( ( kbuoy(:,:,k) == k ) .and. & 0417 ( kbuoy(:,:,k-1) == kxp ) .and. & 0418 ( pfull(:,:,k) >= pshalow ) ) 0419 akhsc(:,:,k) = akhsc0 * rhumscl(:,:) * rhumjmp(:,:) 0420 endwhere 0421 ! ------------------- 0422 end do 0423 0424 ! %%%%%%%%%%%%%%%%%%%%%%%% 0425 endif 0426 0427 !======================================================================= 0428 ! --- CONFINE SHALLOW CONVECTION TO ( pshalow <= p <= plcl ) 0429 !======================================================================= 0430 0431 do k = kctopm1,kxm 0432 ! ------------------- 0433 where ( ( pfull(:,:,k) <= pshalow ) .or. & 0434 ( pfull(:,:,k) >= plcl(:,:) ) ) 0435 akhsc(:,:,k+1) = 0.0 0436 endwhere 0437 ! ------------------- 0438 end do 0439 0440 !======================================================================= 0441 end SUBROUTINE SHALLOW_CONV 0442 0443 !####################################################################### 0444 0445 SUBROUTINE MYLCL ( tlparc, qlparc, plparc, phalf, plcl, kbase ) 0446 0447 !======================================================================= 0448 ! ***** COMPUTE LCL ( CLOUD BASE ) 0449 !======================================================================= 0450 !--------------------------------------------------------------------- 0451 ! Arguments (Intent in) 0452 ! tlparc Initial parcel temperature 0453 ! qlparc Initial parcel mixing ratio 0454 ! plparc Initial parcel pressure 0455 ! phalf Pressure at half levels 0456 ! Arguments (Intent out) 0457 ! plcl Pressure at LCL 0458 ! kbase Index of LCL in column 0459 !--------------------------------------------------------------------- 0460 real, intent(in), dimension(:,:) :: tlparc, qlparc, plparc 0461 real, intent(in), dimension(:,:,:) :: phalf 0462 real, intent(out), dimension(:,:) :: plcl 0463 integer, intent(out), dimension(:,:) :: kbase 0464 0465 !--------------------------------------------------------------------- 0466 ! (Intent local) 0467 !--------------------------------------------------------------------- 0468 integer, parameter :: iter_max = 10 0469 real, parameter :: small = 1.0E-2 0470 0471 real, dimension(size(tlparc,1),size(tlparc,2)) :: & 0472 tlcl, tlclo, clclo, esat, esato, desato, xy1, xy2 0473 0474 logical, dimension(size(tlparc,1),size(tlparc,2)) :: & 0475 non_cnvg 0476 0477 integer :: k, kx, n, iter 0478 0479 !======================================================================= 0480 0481 ! --- Index of lowest model level, etc 0482 kx = size( phalf, 3 ) - 1 0483 0484 ! --- Initial guess for temperature at LCL 0485 tlclo = tlparc 0486 0487 ! --- Compute constant factor 0488 clclo = ( 1.0 + d622/qlparc ) / plparc 0489 clclo = kappa*LOG( clclo ) 0490 clclo = EXP( clclo ) 0491 clclo = tlclo * clclo 0492 0493 ! --- Start with all points non-convergent 0494 non_cnvg = .true. 0495 0496 ! $$$$$$$$$$$$$$$$$$$$ 0497 do iter = 1,iter_max 0498 ! $$$$$$$$$$$$$$$$$$$$ 0499 0500 ! --- Compute saturation vapor pressure and derivative 0501 CALL ESCOMP ( tlclo, esato ) 0502 CALL DESCOMP ( tlclo, desato ) 0503 0504 ! --- Compute new guess for temperature at LCL 0505 where (non_cnvg) 0506 xy1 = kappa * clclo * desato 0507 xy2 = omkappa*LOG( esato ) 0508 xy2 = EXP( xy2 ) 0509 tlcl = ( xy1 * tlclo - clclo * esato ) / ( xy1 - xy2 ) 0510 xy2 = abs( tlcl - tlclo ) 0511 end where 0512 0513 ! --- Test for convergence 0514 where (non_cnvg .and. xy2 <= small) 0515 esat = esato 0516 non_cnvg = .false. 0517 endwhere 0518 n = COUNT( non_cnvg ) 0519 0520 if( n .eq. 0 ) go to 1000 0521 0522 ! --- Shift for next iteration 0523 tlclo = tlcl 0524 0525 ! $$$$$$$$$$$$$$$$$$$$ 0526 end do 0527 ! $$$$$$$$$$$$$$$$$$$$ 0528 !del CALL ERROR_MESG ('MYLCL in SHALLOW_CONV_MOD', & 0529 !del 'ITERATION LOOP FOR LCL FAILED', FATAL) 0530 1000 continue 0531 0532 ! --- Compute pressure at LCL 0533 plcl = ( 1.0 + d622 / qlparc ) * esat 0534 0535 ! --- Bound plcl 0536 plcl(:,:) = MAX( plcl(:,:), pshalow ) 0537 plcl(:,:) = MIN( plcl(:,:), plparc(:,:) ) 0538 0539 ! --- Find index of LCL 0540 do k = 2,kx 0541 where ( ( plcl(:,:) >= phalf(:,:,k ) ) .and. & 0542 ( plcl(:,:) <= phalf(:,:,k+1) ) ) 0543 kbase(:,:) = k 0544 end where 0545 end do 0546 0547 !======================================================================= 0548 end SUBROUTINE MYLCL 0549 0550 !####################################################################### 0551 !####################################################################### 0552 end MODULE SHALLOW_CONV_MOD 0553
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine. The LXR team |
|