MITgcm |
|
|
|
|||
File indexing completed on 2018-03-02 18:37:48 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTCb2ea1d2979 Jean*0001 module vert_turb_driver_mod 0002 0003 !----------------------------------------------------------------------- 0004 ! 0005 ! driver for compuing vertical diffusion coefficients 0006 ! 0007 ! choose either: 0008 ! 1) mellor-yamada 2.5 (with tke) 0009 ! 2) non-local K scheme 0010 ! 0011 !----------------------------------------------------------------------- 0012 !---------------- modules --------------------- 0013 0014 use my25_turb_mod, only: my25_turb_init, my25_turb_end, & 0015 my25_turb, tke_surf, tke 0016 0017 use diffusivity_mod, only: diffusivity, molecular_diff 0018 0019 use shallow_conv_mod, only: shallow_conv_init, shallow_conv 0020 0021 !use diag_manager_mod, only: register_diag_field, send_data 0022 0023 !use time_manager_mod, only: time_type, get_time, operator(-) 0024 0025 use gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit 0026 use constants_mod, only: rdgas, rvgas, kappa 0027 0028 !use fms_mod, only: error_mesg, open_namelist_file, file_exist, & 0029 ! check_nml_error, mpp_root_pe, & 0030 ! mpp_pe, close_file, FATAL, stdlog, & 0031 ! write_version_number 0032 0033 implicit none 0034 private 0035 0036 !---------------- interfaces --------------------- 0037 0038 public vert_turb_driver_init, vert_turb_driver_end, vert_turb_driver 0039 0040 !----------------------------------------------------------------------- 0041 !--------------------- version number ---------------------------------- 0042 15388b052e Jean*0043 character(len=128) :: version = '$Id: vert_turb_driver_mod.F90,v 1.2 2017/08/11 20:48:51 jmc Exp $' b2ea1d2979 Jean*0044 character(len=128) :: tag = '$Name: $' 0045 0046 !----------------------------------------------------------------------- 0047 real, parameter :: p00 = 1000.0E2 0048 real, parameter :: p00inv = 1./p00 0049 real, parameter :: d622 = rdgas/rvgas 0050 real, parameter :: d378 = 1.-d622 0051 real, parameter :: d608 = d378/d622 0052 0053 !---------------- private data ------------------- 0054 0055 logical :: do_init = .true. 0056 0057 real :: gust_zi = 1000. ! constant for computed gustiness (meters) 0058 0059 !----------------------------------------------------------------------- 0060 !-------------------- namelist ----------------------------------------- 0061 0062 logical :: do_shallow_conv = .false. 0063 logical :: do_mellor_yamada = .true. 0064 logical :: use_tau = .true. 0065 logical :: do_molecular_diffusion = .false. 0066 0067 character(len=24) :: gust_scheme = 'constant' ! valid schemes are: 0068 ! => 'constant' 0069 ! => 'beljaars' 0070 real :: constant_gust = 1.0 0071 0072 namelist /vert_turb_driver_nml/ do_shallow_conv, do_mellor_yamada, & 0073 gust_scheme, constant_gust, use_tau, & 0074 do_molecular_diffusion 0075 0076 !-------------------- diagnostics fields ------------------------------- 0077 0078 integer :: id_tke, id_lscale, id_lscale_0, id_z_pbl, id_gust, & 0079 id_diff_t, id_diff_m, id_diff_sc, id_z_full, id_z_half,& 0080 id_uwnd, id_vwnd 0081 0082 real :: missing_value = -999. 0083 0084 character(len=9) :: mod_name = 'vert_turb' 0085 0086 !----------------------------------------------------------------------- 0087 0088 contains 0089 0090 !####################################################################### 0091 0092 subroutine vert_turb_driver (is, js, Time, Time_next, dt, frac_land, & 0093 p_half, p_full, z_half, z_full, & 0094 u_star, b_star, rough, & 0095 uu, vv, tt, qq, & 0096 diff_t, diff_m, gust, & 0097 myThid, mask, kbot ) 0098 ! u, v, t, q, um, vm, tm, qm, & 0099 ! udt, vdt, tdt, qdt, & 0100 0101 !----------------------------------------------------------------------- 0102 integer, intent(in) :: is, js 0103 !type(time_type), intent(in) :: Time, Time_next 0104 real, intent(in) :: Time, Time_next 0105 real, intent(in) :: dt 0106 real, intent(in), dimension(:,:) :: frac_land, & 0107 u_star, b_star, rough 0108 real, intent(in), dimension(:,:,:) :: p_half, p_full, & 0109 z_half, z_full, & 0110 uu, vv, tt, qq 0111 ! u, v, t, q, um, vm, tm, qm, & 0112 ! udt, vdt, tdt, qdt 0113 real, intent(out), dimension(:,:,:) :: diff_t, diff_m 0114 real, intent(out), dimension(:,:) :: gust 0115 integer, intent(in) :: myThid 0116 real, intent(in),optional, dimension(:,:,:) :: mask 0117 integer, intent(in),optional, dimension(:,:) :: kbot 0118 !----------------------------------------------------------------------- 0119 real , dimension(size(tt,1),size(tt,2),size(tt,3)) :: ape, thv 0120 !logical, dimension(size(tt,1),size(tt,2),size(tt,3)+1) :: lmask 0121 !real , dimension(size(tt,1),size(tt,2),size(tt,3)+1) :: diag3 0122 real , dimension(size(tt,1),size(tt,2),size(tt,3)+1) :: el 0123 real , dimension(size(tt,1),size(tt,2)) :: el0, z_pbl 0124 real , dimension(size(diff_t,1),size(diff_t,2), & 0125 size(diff_t,3)) :: diff_sc 0126 !real , dimension(size(t,1),size(t,2),size(t,3)) :: tt, qq, uu, vv 0127 real :: dt_tke 0128 integer :: ie, je, nlev, sec, day 0129 logical :: used 0130 !----------------------------------------------------------------------- 0131 !----------------------- vertical turbulence --------------------------- 0132 !----------------------------------------------------------------------- 0133 0134 if (do_init) CALL PRINT_ERROR( & 0135 'vert_turb_driver in vert_turb_driver_mod'// & 0136 'initialization has not been called', myThid ) 0137 ! if (do_init) call error_mesg & 0138 ! ('vert_turb_driver in vert_turb_driver_mod', & 0139 ! 'initialization has not been called', FATAL) 0140 0141 nlev = size(p_full,3) 0142 ie = is + size(p_full,1) - 1 0143 je = js + size(p_full,2) - 1 0144 0145 !----------------------------------------------------------------------- 0146 !---- set up state variable used by this module ---- 0147 0148 ! if (use_tau) then 0149 ! !-- variables at time tau 0150 ! uu = u 0151 ! vv = v 0152 ! tt = t 0153 ! qq = q 0154 ! else 0155 ! !-- variables at time tau+1 0156 ! uu = um + dt*udt 0157 ! vv = vm + dt*vdt 0158 ! tt = tm + dt*tdt 0159 ! qq = qm + dt*qdt 0160 ! endif 0161 0162 !----------------------------------------------------------------------- 0163 0164 !--------------------------- 0165 if (do_mellor_yamada) then 0166 !--------------------------- 0167 0168 ! ----- time step for prognostic tke calculation ----- 0169 ! call get_time (Time_next-Time, sec, day) 0170 ! dt_tke = real(sec+day*86400) 0171 dt_tke = dt 0172 0173 ! ----- virtual temp ---------- 0174 ape(:,:,:)=(p_full(:,:,:)*p00inv)**(-kappa) 0175 thv(:,:,:)=tt(:,:,:)*(qq(:,:,:)*d608+1.0)*ape(:,:,:) 0176 if (present(mask)) where (mask < 0.5) thv = 200. 0177 0178 ! --------------------- update tke----------------------------------- 0179 ! ---- compute surface tke -------- 0180 ! ---- compute tke, master length scale (el0), ------------- 0181 ! ---- length scale (el), and vert mix coeffs (diff_t,diff_m) ---- 0182 0183 call tke_surf (u_star, tke(is:ie,js:je,:), kbot=kbot) 0184 0185 if ( id_z_pbl > 0 ) then 0186 !------ compute pbl depth from k_profile if diagnostic needed ----- 0187 call my25_turb (dt_tke, frac_land, p_half, p_full, thv, uu, vv, & 0188 z_half, z_full, rough, tke(is:ie,js:je,:), & 0189 el0, el, diff_m, diff_t, & 0190 mask=mask, kbot=kbot, & 0191 ustar=u_star,bstar=b_star,h=z_pbl) 0192 else 0193 call my25_turb (dt_tke, frac_land, p_half, p_full, thv, uu, vv, & 0194 z_half, z_full, rough, tke(is:ie,js:je,:), & 0195 el0, el, diff_m, diff_t, & 0196 mask=mask, kbot=kbot) 0197 end if 0198 0199 !--------------------------- 0200 else 0201 !-------------------------------------------------------------------- 0202 !----------- compute molecular diffusion, if desired --------------- 0203 0204 if (do_molecular_diffusion) then 0205 call molecular_diff (tt, p_half, diff_m, diff_t, myThid ) 0206 else 0207 diff_m = 0.0 0208 diff_t = 0.0 0209 endif 0210 0211 !--------------------------- 0212 !------------------- non-local K scheme -------------- 0213 0214 call diffusivity ( tt, qq, uu, vv, p_full, p_half, z_full, z_half, & 0215 u_star, b_star, z_pbl, diff_m, diff_t, myThid, & 0216 kbot = kbot) 0217 0218 !--------------------------- 0219 endif 0220 !----------------------------------------------------------------------- 0221 !------------------ shallow convection ???? ---------------------------- 0222 0223 if (do_shallow_conv) then 0224 call shallow_conv (tt, qq, p_full, p_half, diff_sc, kbot) 0225 diff_t = diff_t + diff_sc 0226 endif 0227 0228 !----------------------------------------------------------------------- 0229 !------------- define gustiness ------------ 0230 0231 if ( trim(gust_scheme) == 'constant' ) then 0232 gust = constant_gust 0233 else if ( trim(gust_scheme) == 'beljaars' ) then 0234 ! --- from Beljaars (1994) and Beljaars and Viterbo (1999) --- 0235 where (b_star > 0.) 0236 gust = (u_star*b_star*gust_zi)**(1./3.) 0237 elsewhere 0238 gust = 0. 0239 endwhere 0240 endif 0241 0242 !----------------------------------------------------------------------- 0243 !------------------------ diagnostics section -------------------------- 0244 0245 if (do_mellor_yamada) then 0246 0247 ! --- set up local mask for fields with surface data --- 0248 if ( present(mask) ) then 0249 ! lmask(:,:,1) = .true. 0250 ! lmask(:,:,2:nlev+1) = mask(:,:,1:nlev) > 0.5 0251 else 0252 ! lmask = .true. 0253 endif 0254 0255 !------- tke -------------------------------- 0256 if ( id_tke > 0 ) then 0257 ! used = send_data ( id_tke, tke(is:ie,js:je,:), Time_next, & 0258 ! is, js, 1, mask=lmask ) 0259 endif 0260 0261 !------- length scale (at half levels) ------ 0262 if ( id_lscale > 0 ) then 0263 ! used = send_data ( id_lscale, el, Time_next, is, js, 1, & 0264 ! mask=lmask ) 0265 endif 0266 0267 !------- master length scale ------- 0268 if ( id_lscale_0 > 0 ) then 0269 ! used = send_data ( id_lscale_0, el0, Time_next, is, js ) 0270 endif 0271 0272 end if 0273 0274 !------- boundary layer depth ------- 0275 if ( id_z_pbl > 0 ) then 0276 ! used = send_data ( id_z_pbl, z_pbl, Time_next, is, js ) 0277 endif 0278 0279 !------- gustiness ------- 0280 if ( id_gust > 0 ) then 0281 ! used = send_data ( id_gust, gust, Time_next, is, js ) 0282 endif 0283 0284 !------- output diffusion coefficients --------- 0285 0286 if ( id_diff_t > 0 .or. id_diff_m > 0 .or. id_diff_sc > 0 ) then 0287 ! --- set up local mask for fields without surface data --- 0288 if (present(mask)) then 0289 ! lmask(:,:,1:nlev) = mask(:,:,1:nlev) > 0.5 0290 ! lmask(:,:,nlev+1) = .false. 0291 else 0292 ! lmask(:,:,1:nlev) = .true. 0293 ! lmask(:,:,nlev+1) = .false. 0294 endif 0295 ! -- dummy data at surface -- 0296 ! diag3(:,:,nlev+1)=0.0 0297 endif 0298 0299 !------- diffusion coefficient for heat/moisture ------- 0300 if ( id_diff_t > 0 ) then 0301 ! diag3(:,:,1:nlev) = diff_t(:,:,1:nlev) 0302 ! used = send_data ( id_diff_t, diag3, Time_next, is, js, 1, mask=lmask ) 0303 endif 0304 0305 !------- diffusion coefficient for momentum ------- 0306 if ( id_diff_m > 0 ) then 0307 ! diag3(:,:,1:nlev) = diff_m(:,:,1:nlev) 0308 ! used = send_data ( id_diff_m, diag3, Time_next, is, js, 1, mask=lmask ) 0309 endif 0310 0311 !------- diffusion coefficient for shallow conv ------- 0312 if (do_shallow_conv) then 0313 if ( id_diff_sc > 0 ) then 0314 ! diag3(:,:,1:nlev) = diff_sc(:,:,1:nlev) 0315 ! used = send_data ( id_diff_sc, diag3, Time_next, is, js, 1, mask=lmask) 0316 endif 0317 endif 0318 0319 !--- geopotential height relative to the surface on full and half levels ---- 0320 0321 if ( id_z_half > 0 ) then 0322 !--- set up local mask for fields with surface data --- 0323 if ( present(mask) ) then 0324 ! lmask(:,:,1) = .true. 0325 ! lmask(:,:,2:nlev+1) = mask(:,:,1:nlev) > 0.5 0326 else 0327 ! lmask = .true. 0328 endif 0329 ! used = send_data ( id_z_half, z_half, Time_next, is, js, 1, mask=lmask ) 0330 endif 0331 0332 if ( id_z_full > 0 ) then 0333 ! used = send_data ( id_z_full, z_full, Time_next, is, js, 1, rmask=mask) 0334 endif 0335 0336 !--- zonal and meridional wind on mass grid ------- 0337 0338 if ( id_uwnd > 0 ) then 0339 ! used = send_data ( id_uwnd, uu, Time_next, is, js, 1, rmask=mask) 0340 endif 0341 0342 if ( id_vwnd > 0 ) then 0343 ! used = send_data ( id_vwnd, vv, Time_next, is, js, 1, rmask=mask) 0344 endif 0345 0346 !----------------------------------------------------------------------- 0347 0348 end subroutine vert_turb_driver 0349 0350 !####################################################################### 0351 0352 subroutine vert_turb_driver_init (id, jd, kd, axes, Time, myThid) 0353 0354 !----------------------------------------------------------------------- 0355 integer, intent(in) :: id, jd, kd, axes(4) 0356 ! type(time_type), intent(in) :: Time 0357 real, intent(in) :: Time 0358 integer, intent(in) :: myThid 0359 !----------------------------------------------------------------------- 0360 integer, dimension(3) :: full = (/1,2,3/), half = (/1,2,4/) 0361 integer :: ierr, unit, io 0362 0363 integer :: iUnit 0364 CHARACTER*(gcm_LEN_MBUF) :: msgBuf 0365 0366 if (.not.do_init) CALL PRINT_ERROR( & 0367 'vert_turb_driver_init in vert_turb_driver_mod'// & 0368 'attempting to call initialization twice', myThid ) 0369 ! if (.not.do_init) & 0370 ! call error_mesg & 0371 ! ('vert_turb_driver_init in vert_turb_driver_mod', & 0372 ! 'attempting to call initialization twice', FATAL) 0373 0374 !----------------------------------------------------------------------- 0375 !--------------- read namelist ------------------ 0376 0377 ! _BARRIER 0378 ! _BEGIN_MASTER(myThid) 0379 CALL BARRIER(myThid) 0380 IF ( myThid.EQ.1 ) THEN 0381 0382 WRITE(msgBuf,'(A)') 'VERT_TURB_DRIVER_INIT: opening data.atm_gray' 0383 CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid ) 0384 CALL OPEN_COPY_DATA_FILE( & 0385 'data.atm_gray', 'VERT_TURB_DRIVER_INIT', & 0386 iUnit, & 0387 myThid ) 0388 ! Read parameters from open data file 0389 READ(UNIT=iUnit,NML=vert_turb_driver_nml) 0390 WRITE(msgBuf,'(A)') & 0391 'VERT_TURB_DRIVER_INIT: finished reading data.atm_gray' 0392 CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid ) 0393 ! Close the open data file 15388b052e Jean*0394 #ifdef SINGLE_DISK_IO b2ea1d2979 Jean*0395 CLOSE(iUnit) 15388b052e Jean*0396 #else 0397 CLOSE(iUnit,STATUS='DELETE') 0398 #endif /* SINGLE_DISK_IO */ b2ea1d2979 Jean*0399 0400 ! if (file_exist('input.nml')) then 0401 ! unit = open_namelist_file () 0402 ! ierr=1; do while (ierr /= 0) 0403 ! read (unit, nml=vert_turb_driver_nml, iostat=io, end=10) 0404 ! ierr = check_nml_error (io, 'vert_turb_driver_nml') 0405 ! enddo 0406 ! 10 call close_file (unit) 0407 ! endif 0408 0409 !---------- output namelist -------------------------------------------- 0410 0411 ! call write_version_number(version, tag) 0412 ! if ( mpp_pe() == mpp_root_pe() ) write (stdlog(), nml=vert_turb_driver_nml) 0413 0414 ! --- check namelist option --- 0415 if ( trim(gust_scheme) /= 'constant' .and. & 0416 trim(gust_scheme) /= 'beljaars' ) CALL PRINT_ERROR( & 0417 'vert_turb_driver_mod'//'invalid value for namelist'// & 0418 ' variable GUST_SCHEME', myThid ) 0419 ! trim(gust_scheme) /= 'beljaars' ) call error_mesg & 0420 ! ('vert_turb_driver_mod', 'invalid value for namelist & 0421 ! &variable GUST_SCHEME', FATAL) 0422 0423 if (do_molecular_diffusion .and. do_mellor_yamada) CALL PRINT_ERROR( & 0424 'vert_turb_driver_mod'//' cannot activate'// & 0425 ' molecular diffusion with mellor_yamada', myThid ) 0426 ! if (do_molecular_diffusion .and. do_mellor_yamada) & 0427 ! call error_mesg ( 'vert_turb_driver_mod', 'cannot activate & 0428 ! &molecular diffusion with mellor_yamada', FATAL) 0429 0430 !----------------------------------------------------------------------- 0431 0432 if (do_mellor_yamada) call my25_turb_init (id, jd, kd) 0433 0434 if (do_shallow_conv) call shallow_conv_init (kd) 0435 0436 !----------------------------------------------------------------------- 0437 !----- initialize diagnostic fields ----- 0438 0439 id_tke = 0 0440 id_lscale = 0 0441 id_lscale_0 = 0 0442 id_z_pbl = 0 0443 id_gust = 0 0444 id_diff_t = 0 0445 id_diff_m = 0 0446 id_diff_sc = 0 0447 id_z_full = 0 0448 id_z_half = 0 0449 id_uwnd = 0 0450 id_vwnd = 0 0451 ! id_uwnd = register_diag_field ( mod_name, 'uwnd', axes(full), Time, & 0452 ! 'zonal wind on mass grid', 'meters/second' , & 0453 !missing_value=missing_value ) 0454 0455 ! id_vwnd = register_diag_field ( mod_name, 'vwnd', axes(full), Time, & 0456 ! 'meridional wind on mass grid', 'meters/second' , & 0457 !missing_value=missing_value ) 0458 0459 ! id_z_full = & 0460 ! register_diag_field ( mod_name, 'z_full', axes(full), Time, & 0461 ! 'geopotential height relative to surface at full levels', & 0462 ! 'meters' , missing_value=missing_value ) 0463 0464 ! id_z_half = & 0465 ! register_diag_field ( mod_name, 'z_half', axes(half), Time, & 0466 ! 'geopotential height relative to surface at half levels', & 0467 ! 'meters' , missing_value=missing_value ) 0468 0469 if (do_mellor_yamada) then 0470 0471 ! id_tke = & 0472 ! register_diag_field ( mod_name, 'tke', axes(half), Time, & 0473 ! 'turbulent kinetic energy', 'm2/s2' , & 0474 ! missing_value=missing_value ) 0475 0476 ! id_lscale = & 0477 ! register_diag_field ( mod_name, 'lscale', axes(half), Time, & 0478 ! 'turbulent length scale', 'm' , & 0479 ! missing_value=missing_value ) 0480 0481 ! id_lscale_0 = & 0482 ! register_diag_field ( mod_name, 'lscale_0', axes(1:2), Time, & 0483 ! 'master length scale', 'm' ) 0484 endif 0485 0486 ! id_z_pbl = & 0487 ! register_diag_field ( mod_name, 'z_pbl', axes(1:2), Time, & 0488 ! 'depth of planetary boundary layer', 'm' ) 0489 0490 ! id_gust = & 0491 ! register_diag_field ( mod_name, 'gust', axes(1:2), Time, & 0492 ! 'wind gustiness in surface layer', 'm/s' ) 0493 0494 ! id_diff_t = & 0495 ! register_diag_field ( mod_name, 'diff_t', axes(half), Time, & 0496 ! 'vert diff coeff for temp', 'm2/s' , & 0497 ! missing_value=missing_value ) 0498 0499 ! id_diff_m = & 0500 ! register_diag_field ( mod_name, 'diff_m', axes(half), Time, & 0501 ! 'vert diff coeff for momentum', 'm2/s' , & 0502 ! missing_value=missing_value ) 0503 0504 if (do_shallow_conv) then 0505 0506 ! id_diff_sc = & 0507 ! register_diag_field ( mod_name, 'diff_sc', axes(half), Time, & 0508 ! 'vert diff coeff for shallow conv', 'm2/s' , & 0509 ! missing_value=missing_value ) 0510 endif 0511 0512 !----------------------------------------------------------------------- 0513 0514 do_init=.false. 0515 0516 ENDIF 0517 CALL BARRIER(myThid) 0518 0519 !----------------------------------------------------------------------- 0520 0521 end subroutine vert_turb_driver_init 0522 0523 !####################################################################### 0524 0525 subroutine vert_turb_driver_end(myThid) 0526 0527 integer, intent(in) :: myThid 0528 0529 !----------------------------------------------------------------------- 0530 if (do_mellor_yamada) call my25_turb_end 0531 !----------------------------------------------------------------------- 0532 0533 end subroutine vert_turb_driver_end 0534 0535 !####################################################################### 0536 0537 end module vert_turb_driver_mod
| [ 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 |
|