Index: LMDZ6/trunk/DefLists/field_def_lmdz.xml
===================================================================
--- LMDZ6/trunk/DefLists/field_def_lmdz.xml (revision 5036)
+++ LMDZ6/trunk/DefLists/field_def_lmdz.xml (revision 5039)
@@ -636,4 +636,5 @@
+
Index: LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90 (revision 5036)
+++ LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90 (revision 5039)
@@ -7,5 +7,5 @@
subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
- tke,eps,Km_out,Kh_out)
+ tke,eps,tke_shear,tke_buoy,tke_trans,Km_out,Kh_out)
!========================================================================
@@ -79,4 +79,7 @@
REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: eps ! output: TKE dissipation rate at interface between layers (m2/s3)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_shear! output: TKE shear production rate (m2/s3)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_buoy ! output: TKE buoyancy production rate (m2/s3)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_trans! output: TKE transport (diffusion) term (m2/s3)
REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers (m2/s)
REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers (m2/s)
@@ -261,4 +264,7 @@
shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
+ tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay)
+ tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) &
+ *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
ENDDO
ENDDO
@@ -278,5 +284,8 @@
qq=max(0.,qq)
tke(igrid,ilay)=0.5*(qq**2)
- eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
+ eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
+ tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay)
+ tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) &
+ *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
ENDDO
ENDDO
@@ -293,5 +302,8 @@
qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.) &
*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
- /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
+ /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
+ tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
+ tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
+ *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
tke(igrid,ilay)=0.5*(qq**2)
eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
@@ -308,4 +320,7 @@
eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay))
qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+ tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
+ tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
+ *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
IF (Ri(igrid,ilay) .LT. 0.) THEN
netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
@@ -327,4 +342,7 @@
DO igrid=1,ngrid
qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
+ tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay)
+ tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) &
+ *(Ri(igrid,ilay) / Prandtl(igrid,ilay))
qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
+qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
@@ -349,4 +367,6 @@
tke(igrid,nlay+1)=0.
eps(igrid,nlay+1)=0.
+ tke_shear(igrid,nlay+1)=0.
+ tke_buoy(igrid,nlay+1)=0.
END DO
@@ -359,4 +379,6 @@
tke(igrid,1)=ctkes*(ustar**2)
eps(igrid,1)=0. ! arbitrary as TKE is not properly defined at the surface
+ tke_shear(igrid,1)=0.
+ tke_buoy(igrid,1)=0.
END DO
@@ -364,6 +386,7 @@
! vertical diffusion of TKE
!==========================
+tke_trans(:,:)=0.
IF (atke_ok_vdiff) THEN
- CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
+ CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans)
ENDIF
@@ -387,5 +410,5 @@
!===============================================================================================
-subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
+subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans)
! routine that computes the vertical diffusion of TKE by the turbulence
@@ -408,5 +431,5 @@
REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers
-
+REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke_trans ! turbulent kinetic energy transport term (m2/s3)
@@ -480,4 +503,5 @@
! update TKE
tke(:,:)=tke(:,:)+dtke(:,:)
+tke_trans(:,:)=dtke(:,:)/dtime
Index: LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90 (revision 5036)
+++ LMDZ6/trunk/libf/phylmd/lmdz_call_atke.F90 (revision 5039)
@@ -8,5 +8,5 @@
contains
-subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
+subroutine call_atke(dtime,ngrid,nlay,nsrf,ni,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
wind_u,wind_v,temp,qvap,play,pinterf, &
tke,eps,Km_out,Kh_out)
@@ -16,4 +16,5 @@
USE lmdz_atke_turbulence_ini, ONLY : iflag_num_atke, rg, rd
+USE phys_local_var_mod, ONLY: tke_shear, tke_buoy, tke_trans
implicit none
@@ -26,4 +27,6 @@
INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
INTEGER, INTENT(IN) :: nlay ! number of vertical index
+INTEGER, INTENT(IN) :: nsrf ! surface tile index
+INTEGER, DIMENSION(ngrid), INTENT(IN) :: ni ! array of indices to move from knon to klon arrays
@@ -50,12 +53,13 @@
+REAL, DIMENSION(ngrid,nlay+1) :: tke_shear_term,tke_buoy_term,tke_trans_term
REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict
REAL, DIMENSION(ngrid) :: wind1
-INTEGER i
+INTEGER i,j,k
call atke_compute_km_kh(ngrid,nlay,dtime,&
wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv,&
- tke,eps,Km_out,Kh_out)
+ tke,eps,tke_shear_term,tke_buoy_term,tke_trans_term,Km_out,Kh_out)
@@ -76,10 +80,19 @@
call atke_compute_km_kh(ngrid,nlay,dtime,&
wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,cdrag_uv, &
- tke,eps,Km_out,Kh_out)
+ tke,eps,tke_shear_term,tke_buoy_term,tke_trans_term,Km_out,Kh_out)
end if
+! Diagnostics of tke loss/source terms
+ DO k=1,nlay+1
+ DO i=1,ngrid
+ j=ni(i)
+ tke_shear(j,k,nsrf)=tke_shear_term(i,k)
+ tke_buoy(j,k,nsrf)=tke_buoy_term(i,k)
+ tke_trans(j,k,nsrf)=tke_trans_term(i,k)
+ ENDDO
+ ENDDO
Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 5036)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 5039)
@@ -1996,5 +1996,5 @@
IF (iflag_pbl>=50) THEN
- CALL call_atke(dtime,knon,klev,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &
+ CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &
yu(1:knon,:),yv(1:knon,:),yt(1:knon,:),yq(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &
ytke(1:knon,:),yeps(1:knon,:), ycoefm(1:knon,:), ycoefh(1:knon,:))
@@ -2041,5 +2041,5 @@
IF (iflag_pbl>=50) THEN
- CALL call_atke(dtime,knon,klev,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon), &
+ CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon), &
yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &
ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))
@@ -2081,5 +2081,5 @@
IF (iflag_pbl>=50) THEN
- CALL call_atke(dtime,knon,klev,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &
+ CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &
yu_w(1:knon,:),yv_w(1:knon,:),yt_w(1:knon,:),yq_w(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:), &
ytke_w(1:knon,:),yeps_w(1:knon,:),ycoefm_w(1:knon,:),ycoefh_w(1:knon,:))
Index: LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 5036)
+++ LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 5039)
@@ -26,4 +26,6 @@
REAL, SAVE, ALLOCATABLE :: pbl_eps(:,:,:)
!$OMP THREADPRIVATE(pbl_eps)
+ REAL, SAVE, ALLOCATABLE :: tke_shear(:,:,:), tke_buoy(:,:,:), tke_trans(:,:,:)
+ !$OMP THREADPRIVATE(tke_shear(:,:,:),tke_buoy(:,:,:),tke_trans(:,:,:))
REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
!$OMP THREADPRIVATE(tr_seri)
@@ -694,5 +696,7 @@
ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
+ ALLOCATE(tke_shear(klon,klev+1,nbsrf), tke_buoy(klon,klev+1,nbsrf), tke_trans(klon,klev+1,nbsrf))
pbl_eps(:,:,:)=0.
+ tke_shear(:,:,:)=0.; tke_buoy(:,:,:)=0.; tke_trans(:,:,:)=0.
l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
ALLOCATE(rhcl(klon,klev))
@@ -1055,4 +1059,5 @@
DEALLOCATE(u_seri,v_seri)
DEALLOCATE(l_mixmin,l_mix,wprime)
+ DEALLOCATE(tke_shear,tke_buoy,tke_trans)
DEALLOCATE(pbl_eps)
DEALLOCATE(rhcl)
Index: LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 5036)
+++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 5039)
@@ -1114,6 +1114,13 @@
TYPE(ctrl_out), SAVE :: o_tke = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
'tke ', 'TKE', 'm2/s2', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_tke_shear = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'tke_shear ', 'TKE shear term', 'm2/s3', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_tke_buoy = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'tke_buoy ', 'TKE buoyancy term', 'm2/s3', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_tke_trans = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'tke_trans ', 'TKE transport term', 'm2/s3', (/ ('', i=1, 10) /))
TYPE(ctrl_out), SAVE :: o_tke_dissip = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
- 'tke_dissip ', 'TKE DISSIPATION', 'm2/s3', (/ ('', i=1, 10) /))
+ 'tke_dissip ', 'TKE dissipation term', 'm2/s3', (/ ('', i=1, 10) /))
+
TYPE(ctrl_out), SAVE :: o_tke_max = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
'tke_max', 'TKE max', 'm2/s2', &
Index: LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 5036)
+++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 5039)
@@ -146,5 +146,5 @@
o_dqsphy, o_dqsphy2d, o_dqbsphy, o_dqbsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
o_ages_srf, o_snow_srf, o_alb1, o_alb2, o_tke, o_tke_dissip, &
- o_tke_max, o_kz, o_kz_max, o_clwcon, &
+ o_tke_max, o_kz, o_kz_max, o_clwcon, o_tke_shear, o_tke_buoy, o_tke_trans, &
o_dtdyn, o_dqdyn, o_dqdyn2d, o_dqldyn, o_dqldyn2d, &
o_dqsdyn, o_dqsdyn2d, o_dqbsdyn, o_dqbsdyn2d, o_dudyn, o_dvdyn, &
@@ -314,5 +314,5 @@
zn2mout, t2m_min_mon, t2m_max_mon, evap, &
snowerosion, zxustartlic, zxrhoslic, zxqsaltlic, &
- l_mixmin,l_mix, pbl_eps, &
+ l_mixmin,l_mix, pbl_eps, tke_shear, tke_buoy, tke_trans, &
zu10m, zv10m, zq2m, zustar, zxqsurf, &
rain_lsc, rain_num, snow_lsc, bils, sens, fder, &
@@ -1313,5 +1313,43 @@
ENDIF
- CALL histwrite_phy(o_tke_dissip, zx_tmp_fi3d)
+ CALL histwrite_phy(o_tke_dissip, zx_tmp_fi3d)
+
+ zx_tmp_fi3d=0.
+ IF (vars_defined) THEN
+ DO nsrf=1,nbsrf
+ DO k=1,klev
+ zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) &
+ +pctsrf(:,nsrf)*tke_shear(:,k,nsrf)
+ ENDDO
+ ENDDO
+ ENDIF
+
+ CALL histwrite_phy(o_tke_shear, zx_tmp_fi3d)
+
+ zx_tmp_fi3d=0.
+ IF (vars_defined) THEN
+ DO nsrf=1,nbsrf
+ DO k=1,klev
+ zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) &
+ +pctsrf(:,nsrf)*tke_buoy(:,k,nsrf)
+ ENDDO
+ ENDDO
+ ENDIF
+
+ CALL histwrite_phy(o_tke_buoy, zx_tmp_fi3d)
+
+
+ zx_tmp_fi3d=0.
+ IF (vars_defined) THEN
+ DO nsrf=1,nbsrf
+ DO k=1,klev
+ zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) &
+ +pctsrf(:,nsrf)*tke_trans(:,k,nsrf)
+ ENDDO
+ ENDDO
+ ENDIF
+
+ CALL histwrite_phy(o_tke_trans, zx_tmp_fi3d)
+
ENDIF
Index: LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90 (revision 5036)
+++ LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90 (revision 5039)
@@ -28,4 +28,6 @@
REAL, SAVE, ALLOCATABLE :: pbl_eps(:,:,:)
!$OMP THREADPRIVATE(pbl_eps)
+ REAL, SAVE, ALLOCATABLE :: tke_shear(:,:,:), tke_buoy(:,:,:), tke_trans(:,:,:)
+ !$OMP THREADPRIVATE(tke_shear(:,:,:),tke_buoy(:,:,:),tke_trans(:,:,:))
REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
!$OMP THREADPRIVATE(tr_seri)
@@ -815,5 +817,7 @@
ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
+ ALLOCATE(tke_shear(klon,klev+1,nbsrf), tke_buoy(klon,klev+1,nbsrf), tke_trans(klon,klev+1,nbsrf))
pbl_eps(:,:,:)=0.
+ tke_shear(:,:,:)=0.; tke_buoy(:,:,:)=0.; tke_trans(:,:,:)=0.
l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
ALLOCATE(rhcl(klon,klev))
@@ -1253,4 +1257,5 @@
DEALLOCATE(u_seri,v_seri)
DEALLOCATE(l_mixmin,l_mix,wprime)
+ DEALLOCATE(tke_shear,tke_buoy,tke_trans)
DEALLOCATE(pbl_eps)
DEALLOCATE(rhcl)
Index: LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90 (revision 5036)
+++ LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90 (revision 5039)
@@ -1114,6 +1114,13 @@
TYPE(ctrl_out), SAVE :: o_tke = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
'tke ', 'TKE', 'm2/s2', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_tke_shear = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'tke_shear ', 'TKE shear term', 'm2/s3', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_tke_buoy = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'tke_buoy ', 'TKE buoyancy term', 'm2/s3', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_tke_trans = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+ 'tke_trans ', 'TKE transport term', 'm2/s3', (/ ('', i=1, 10) /))
TYPE(ctrl_out), SAVE :: o_tke_dissip = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
- 'tke_dissip ', 'TKE DISSIPATION', 'm2/s3', (/ ('', i=1, 10) /))
+ 'tke_dissip ', 'TKE dissipation term', 'm2/s3', (/ ('', i=1, 10) /))
+
TYPE(ctrl_out), SAVE :: o_tke_max = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
'tke_max', 'TKE max', 'm2/s2', &