Index: LMDZ6/trunk/libf/phylmdiso/albsno.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/albsno.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/albsno.f90 (revision 5986) @@ -0,0 +1,62 @@ +! +! $Header$ +! +SUBROUTINE albsno(klon, knon, dtime, agesno, alb_neig_grid, precip_snow) + + USE clesphys_mod_h + IMPLICIT NONE + + +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: klon, knon + REAL, INTENT(IN) :: dtime + REAL, DIMENSION(klon), INTENT(IN) :: precip_snow + +! In/Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(INOUT) :: agesno + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(OUT) :: alb_neig_grid + +! Local variables +!**************************************************************************************** + INTEGER :: i, nv + INTEGER, PARAMETER :: nvm = 8 + REAL :: as + REAL, DIMENSION(klon,nvm) :: veget + REAL, DIMENSION(nvm),SAVE :: init, decay + !$OMP THREADPRIVATE(init, decay) + + DATA init /0.55, 0.14, 0.18, 0.29, 0.15, 0.15, 0.14, 0./ + DATA decay/0.30, 0.67, 0.63, 0.45, 0.40, 0.14, 0.06, 1./ +!**************************************************************************************** + + if (albsno0>=0.) then + init(:)=albsno0 + decay(:)=0. + endif + + veget = 0. + veget(:,1) = 1. ! desert partout + DO i = 1, knon + alb_neig_grid(i) = 0.0 + ENDDO + DO nv = 1, nvm + DO i = 1, knon + as = init(nv)+decay(nv)*EXP(-agesno(i)/5.) + alb_neig_grid(i) = alb_neig_grid(i) + veget(i,nv)*as + ENDDO + ENDDO + + +! modilation en fonction de l'age de la neige + DO i = 1, knon + agesno(i) = (agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)& + & * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/0.3) + agesno(i) = MAX(agesno(i),0.0) + ENDDO + +END SUBROUTINE albsno Index: LMDZ6/trunk/libf/phylmdiso/calbeta.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/calbeta.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/calbeta.f90 (revision 5986) @@ -0,0 +1,113 @@ +! +! $Header$ +! +SUBROUTINE calbeta(dtime,indice,knon,snow,qsol, & + vbeta,vcal,vdif) + +USE flux_arp_mod_h + USE dimphy + USE indice_sol_mod + + IMPLICIT none + + +!====================================================================== +! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM au LMD) +! date: 19940414 +!====================================================================== +! +! Calculer quelques parametres pour appliquer la couche limite +! ------------------------------------------------------------ +! Variables d'entrees +!**************************************************************************************** + REAL, INTENT(IN) :: dtime + INTEGER, INTENT(IN) :: indice + INTEGER, INTENT(IN) :: knon + REAL, DIMENSION(klon), INTENT(IN) :: snow + REAL, DIMENSION(klon), INTENT(IN) :: qsol + + +! Variables de sorties +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(OUT) :: vbeta + REAL, DIMENSION(klon), INTENT(OUT) :: vcal + REAL, DIMENSION(klon), INTENT(OUT) :: vdif + +! Variables locales +!**************************************************************************************** + REAL, PARAMETER :: tau_gl=86400.0*5.0 ! temps de relaxation pour la glace de mer +!cc PARAMETER (tau_gl=86400.0*30.0) + REAL, PARAMETER :: mx_eau_sol=150.0 + REAL, PARAMETER :: calsol=1.0/(2.5578E+06*0.15) + REAL, PARAMETER :: calsno=1.0/(2.3867E+06*0.15) + REAL, PARAMETER :: calice=1.0/(5.1444E+06*0.15) + + INTEGER :: i + +!**************************************************************************************** + + vbeta(:) = 0.0 + vcal(:) = 0.0 + vdif(:) = 0.0 + + IF (indice.EQ.is_oce) THEN + DO i = 1, knon + vcal(i) = 0.0 + vbeta(i) = 1.0 + vdif(i) = 0.0 + ENDDO + ENDIF + + IF (indice.EQ.is_sic) THEN + DO i = 1, knon + vcal(i) = calice + IF (snow(i) .GT. 0.0) vcal(i) = calsno + vbeta(i) = 1.0 + vdif(i) = 1.0/tau_gl +! vdif(i) = calice/tau_gl ! c'etait une erreur + ENDDO + ENDIF + + IF (indice.EQ.is_ter) THEN + DO i = 1, knon + vcal(i) = calsol + IF (snow(i) .GT. 0.0) vcal(i) = calsno + vbeta(i) = MIN(2.0*qsol(i)/mx_eau_sol, 1.0) + vdif(i) = 0.0 + ENDDO + ENDIF + + IF (indice.EQ.is_lic) THEN + DO i = 1, knon + vcal(i) = calice + IF (snow(i) .GT. 0.0) vcal(i) = calsno + vbeta(i) = 1.0 + vdif(i) = 0.0 + ENDDO + ENDIF + + ! EV: when beta is prescribed for 1D cases: + IF (knon.EQ.1 .AND. ok_prescr_beta) THEN + DO i = 1, knon + vbeta(i)=betaevap + ENDDO + ENDIF + +END SUBROUTINE calbeta + + + + + + + + + + + + + + + + + Index: LMDZ6/trunk/libf/phylmdiso/calbeta_clim.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/calbeta_clim.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/calbeta_clim.f90 (revision 5986) @@ -0,0 +1,72 @@ +! +! $Header: /home/cvsroot/LMDZ4/libf/phylmd/calbeta.F90,v 1.2 2007/06/22 12:49:51 +! fairhead Exp $ +! + +SUBROUTINE calbeta_clim(klon,time,lat_radian,beta) + + !====================================================================== + ! Auteur(s): A.K. TRAORE + !====================================================================== + + !USE phys_local_var_mod, ONLY : ideal_beta !pour faire la variable dans le + ! physiq.f pour des sorties directes de beta + + USE phys_cal_mod, only: year_len + USE print_control_mod, ONLY: prt_level + + implicit none + integer klon,nt,j,it + real logbeta(klon),pi + real lat(klon),lat_radian(klon) + integer time + real time_radian + real lat_sahel,beta(klon) + real lat_nord,lat_sud + + !============================================== + + pi=2.*asin(1.) + beta=0. + + !calcul des cordonnees + + ! print*,'LATITUDES BETA ',lat_radian + time_radian=(time+15.)*2.*pi / year_len + + if (prt_level >= 1) print *, 'time_radian time', time_radian, time + + lat(:)=180.*lat_radian(:)/pi !lat(:)=lat_radian(:) + + lat_sahel=-5*sin(time_radian)+13 + lat_nord=lat_sahel+25. + lat_sud=lat_sahel-25. + do j=1,klon + !=========== + if (lat(j) < 5. ) then + + logbeta(j)=0.2*(lat(j)-lat_sud)-1.6 + beta(j)=10**(logbeta(j)) + beta(j)=max(beta(j),0.03) + beta(j)=min(beta(j),0.22) + ! print*,'j,lat,lat_radian,beta',j,lat(j),lat_radian(j),beta(j) + !=========== + elseif (lat(j) < 22.) then !lat(j)<22. + + logbeta(j)=-0.25*(lat(j)-lat_sahel)-1.6 + beta(j)=10**(logbeta(j)) + beta(j)=max(beta(j),1.e-2) + beta(j)=min(beta(j),0.22) + ! print*,'j,lat,lat_radian,beta',j,lat(j),lat_radian(j),beta(j) + !=========== + else + logbeta(j)=0.25*(lat(j)-lat_nord)-1. + beta(j)=10**(logbeta(j)) + beta(j)=max(beta(j),1.e-2) + beta(j)=min(beta(j),0.25) + ! print*,'j,lat,lat_radian,beta',j,lat(j),lat_radian(j),beta(j) + endif + !=========== + enddo + +end SUBROUTINE calbeta_clim Index: LMDZ6/trunk/libf/phylmdiso/climb_qbs_mod.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/climb_qbs_mod.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/climb_qbs_mod.f90 (revision 5986) @@ -0,0 +1,379 @@ +MODULE climb_qbs_mod +! +! Module to solve the verctical diffusion of blowing snow; +! + USE dimphy + + IMPLICIT NONE + SAVE + PRIVATE + PUBLIC :: climb_qbs_down, climb_qbs_up + + REAL, DIMENSION(:,:), ALLOCATABLE :: gamaqbs + !$OMP THREADPRIVATE(gamaqbs) + REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_QBS, Dcoef_QBS + !$OMP THREADPRIVATE(Ccoef_QBS, Dcoef_QBS) + REAL, DIMENSION(:), ALLOCATABLE :: Acoef_QBS, Bcoef_QBS + !$OMP THREADPRIVATE(Acoef_QBS, Bcoef_QBS) + REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefqbs + !$OMP THREADPRIVATE(Kcoefqbs) + +CONTAINS +! +!**************************************************************************************** +! + SUBROUTINE climb_qbs_down(knon, coefqbs, paprs, pplay, & + delp, temp, qbs, dtime, & + Ccoef_QBS_out, Dcoef_QBS_out, & + Kcoef_qbs_out, gama_qbs_out, & + Acoef_QBS_out, Bcoef_QBS_out) + +! This routine calculates recursivly the coefficients C and D +! for the quantity X=[QBS] in equation X(k) = C(k) + D(k)*X(k-1), where k is +! the index of the vertical layer. +USE yomcst_mod_h +USE compbl_mod_h +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: knon + REAL, DIMENSION(klon,klev), INTENT(IN) :: coefqbs + REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay + REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs + REAL, DIMENSION(klon,klev), INTENT(IN) :: delp + REAL, DIMENSION(klon,klev), INTENT(IN) :: temp + REAL, DIMENSION(klon,klev), INTENT(IN) :: qbs + REAL, INTENT(IN) :: dtime + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_QBS_out + REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_QBS_out + + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_QBS_out + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_QBS_out + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Kcoef_qbs_out + REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_qbs_out + +! Local variables +!**************************************************************************************** + LOGICAL, SAVE :: first=.TRUE. + !$OMP THREADPRIVATE(first) + REAL, DIMENSION(klon) :: psref + REAL :: delz, pkh + INTEGER :: k, i, ierr +!**************************************************************************************** +! 1) +! Allocation at first time step only +! +!**************************************************************************************** + + IF (first) THEN + first=.FALSE. + ALLOCATE(Ccoef_QBS(klon,klev), STAT=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_QBS, ierr=', ierr + + ALLOCATE(Dcoef_QBS(klon,klev), STAT=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_QBS, ierr=', ierr + + ALLOCATE(Acoef_QBS(klon), Bcoef_QBS(klon), STAT=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_BS and Bcoef_BS, ierr=', ierr + + ALLOCATE(Kcoefqbs(klon,klev), STAT=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in allloc Kcoefqbs, ierr=', ierr + + ALLOCATE(gamaqbs(1:klon,2:klev), STAT=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaqbs, ierr=', ierr + + END IF + +!**************************************************************************************** +! 2) +! Definition of the coeficient K +! +!**************************************************************************************** + Kcoefqbs(:,:) = 0.0 + DO k = 2, klev + DO i = 1, knon + Kcoefqbs(i,k) = & + coefqbs(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) & + *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2 + ENDDO + ENDDO + +!**************************************************************************************** +! 3) +! Calculation of gama for "Q" and "H" +! +!**************************************************************************************** +! surface pressure is used as reference + psref(:) = paprs(:,1) + +! definition of gama + IF (iflag_pbl == 1) THEN + gamaqbs(:,:) = 0.0 + +! conversion de gama + DO k = 2, klev + DO i = 1, knon + delz = RD * (temp(i,k-1)+temp(i,k)) / & + 2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k)) + pkh = (psref(i)/paprs(i,k))**RKAPPA + +! convertie gradient verticale de contenu en neige soufflee en difference de neige soufflee entre centre de couches + gamaqbs(i,k) = gamaqbs(i,k) * delz + ENDDO + ENDDO + + ELSE + gamaqbs(:,:) = 0.0 + ENDIF + + +!**************************************************************************************** +! 4) +! Calculte the coefficients C and D for specific content of blowing snow, qbs +! +!**************************************************************************************** + + CALL calc_coef_qbs(knon, Kcoefqbs(:,:), gamaqbs(:,:), delp(:,:), qbs(:,:), & + Ccoef_QBS(:,:), Dcoef_QBS(:,:), Acoef_QBS, Bcoef_QBS) + + +!**************************************************************************************** +! 5) +! Return the first layer in output variables +! +!**************************************************************************************** + Acoef_QBS_out = Acoef_QBS + Bcoef_QBS_out = Bcoef_QBS + +!**************************************************************************************** +! 6) +! If Pbl is split, return also the other layers in output variables +! +!**************************************************************************************** + IF (mod(iflag_pbl_split,10) .ge.1) THEN + DO k= 1, klev + DO i= 1, klon + Ccoef_QBS_out(i,k) = Ccoef_QBS(i,k) + Dcoef_QBS_out(i,k) = Dcoef_QBS(i,k) + Kcoef_qbs_out(i,k) = Kcoefqbs(i,k) + IF (k.eq.1) THEN + gama_qbs_out(i,k) = 0. + ELSE + gama_qbs_out(i,k) = gamaqbs(i,k) + ENDIF + ENDDO + ENDDO +!!! + ENDIF ! (mod(iflag_pbl_split,2) .ge.1) +!!! + + END SUBROUTINE climb_qbs_down +! +!**************************************************************************************** +! + SUBROUTINE calc_coef_qbs(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef) +! +! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1) +! where X is QQBS, and k the vertical level k=1,klev +USE yomcst_mod_h +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: knon + REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp + REAL, DIMENSION(klon,klev), INTENT(IN) :: X + REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef + +! Local variables +!**************************************************************************************** + INTEGER :: k, i + REAL :: buf + +!**************************************************************************************** +! Niveau au sommet, k=klev +! +!**************************************************************************************** + Ccoef(:,:) = 0.0 + Dcoef(:,:) = 0.0 + + DO i = 1, knon + buf = delp(i,klev) + Kcoef(i,klev) + + Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf + Dcoef(i,klev) = Kcoef(i,klev)/buf + END DO + + +!**************************************************************************************** +! Niveau (klev-1) <= k <= 2 +! +!**************************************************************************************** + + DO k=(klev-1),2,-1 + DO i = 1, knon + buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1)) + Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + & + Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf + Dcoef(i,k) = Kcoef(i,k)/buf + END DO + END DO + +!**************************************************************************************** +! Niveau k=1 +! +!**************************************************************************************** + + DO i = 1, knon + buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2)) + Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf + Bcoef(i) = -1. * RG / buf + END DO + + END SUBROUTINE calc_coef_qbs +! +!**************************************************************************************** +! + SUBROUTINE climb_qbs_up(knon, dtime, qbs_old, & + flx_qbs1, paprs, pplay, & + Acoef_QBS_in, Bcoef_QBS_in, & + Ccoef_QBS_in, Dcoef_QBS_in, & + Kcoef_qbs_in, gama_qbs_in, & + flux_qbs, d_qbs) +! +! This routine calculates the flux and tendency of the specific content of blowing snow qbs +! The quantity qbs is calculated according to +! X(k) = C(k) + D(k)*X(k-1) for X=[qbs], where the coefficients +! C and D are known from before and k is index of the vertical layer. +! +USE yomcst_mod_h +USE compbl_mod_h +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: knon + REAL, INTENT(IN) :: dtime + REAL, DIMENSION(klon,klev), INTENT(IN) :: qbs_old + REAL, DIMENSION(klon), INTENT(IN) :: flx_qbs1 + REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs + REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay + +!!! nrlmd le 02/05/2011 + REAL, DIMENSION(klon), INTENT(IN) :: Acoef_QBS_in, Bcoef_QBS_in + REAL, DIMENSION(klon,klev), INTENT(IN) :: Ccoef_QBS_in, Dcoef_QBS_in + REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef_qbs_in, gama_qbs_in +!!! + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon,klev), INTENT(OUT) :: flux_qbs, d_qbs + +! Local variables +!**************************************************************************************** + LOGICAL, SAVE :: last=.FALSE. +!$OMP THREADPRIVATE(last) + REAL, DIMENSION(klon,klev) :: qbs_new + REAL, DIMENSION(klon) :: psref + INTEGER :: k, i, ierr +!**************************************************************************************** +! 1) +! Definition of some variables + REAL, DIMENSION(klon,klev) :: zairm +! +!**************************************************************************************** + flux_qbs(:,:) = 0.0 + d_qbs(:,:) = 0.0 + + psref(1:knon) = paprs(1:knon,1) + + IF (mod(iflag_pbl_split,10) .ge.1) THEN + DO i = 1, knon + Acoef_QBS(i)=Acoef_QBS_in(i) + Bcoef_QBS(i)=Bcoef_QBS_in(i) + ENDDO + DO k = 1, klev + DO i = 1, knon + Ccoef_QBS(i,k)=Ccoef_QBS_in(i,k) + Dcoef_QBS(i,k)=Dcoef_QBS_in(i,k) + Kcoefqbs(i,k)=Kcoef_qbs_in(i,k) + IF (k.gt.1) THEN + gamaqbs(i,k)=gama_qbs_in(i,k) + ENDIF + ENDDO + ENDDO +!!! + ENDIF ! (mod(iflag_pbl_split,2) .ge.1) +!!! + +!**************************************************************************************** +! 2) +! Calculation of QBS +! +!**************************************************************************************** + +!- First layer + qbs_new(1:knon,1) = Acoef_QBS(1:knon) + Bcoef_QBS(1:knon)*flx_qbs1(1:knon)*dtime +!- All the other layers + DO k = 2, klev + DO i = 1, knon + qbs_new(i,k) = Ccoef_QBS(i,k) + Dcoef_QBS(i,k)*qbs_new(i,k-1) + END DO + END DO +!**************************************************************************************** +! 3) +! Calculation of the flux for QBS +! +!**************************************************************************************** + +!- The flux at first layer, k=1 + flux_qbs(1:knon,1)=flx_qbs1(1:knon) + +!- The flux at all layers above surface + DO k = 2, klev + DO i = 1, knon + flux_qbs(i,k) = (Kcoefqbs(i,k)/RG/dtime) * & + (qbs_new(i,k)-qbs_new(i,k-1)+gamaqbs(i,k)) + END DO + END DO + +!**************************************************************************************** +! 4) +! Calculation of tendency for QBS +! +!**************************************************************************************** + DO k = 1, klev + DO i = 1, knon + d_qbs(i,k) = qbs_new(i,k) - qbs_old(i,k) + zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg + END DO + END DO + +!**************************************************************************************** +! Some deallocations +! +!**************************************************************************************** + IF (last) THEN + DEALLOCATE(Ccoef_QBS, Dcoef_QBS,stat=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_QBS, Dcoef_QBS, ierr=', ierr + DEALLOCATE(Acoef_QBS, Bcoef_QBS,stat=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_QBS, Bcoef_QBS, ierr=', ierr + DEALLOCATE(gamaqbs,stat=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaqbs, ierr=', ierr + DEALLOCATE(Kcoefqbs,stat=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefqbs, ierr=', ierr + END IF + END SUBROUTINE climb_qbs_up +! +!**************************************************************************************** +! +END MODULE climb_qbs_mod + + + + + + Index: LMDZ6/trunk/libf/phylmdiso/climb_wind_mod.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/climb_wind_mod.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/climb_wind_mod.f90 (revision 5986) @@ -0,0 +1,376 @@ +! +MODULE climb_wind_mod +! +! Module to solve the verctical diffusion of the wind components "u" and "v". +! + USE dimphy + + IMPLICIT NONE + + SAVE + PRIVATE + + REAL, DIMENSION(:), ALLOCATABLE :: alf1, alf2 + !$OMP THREADPRIVATE(alf1,alf2) + REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefm + !$OMP THREADPRIVATE(Kcoefm) + REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_U, Dcoef_U + !$OMP THREADPRIVATE(Ccoef_U, Dcoef_U) + REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_V, Dcoef_V + !$OMP THREADPRIVATE(Ccoef_V, Dcoef_V) + REAL, DIMENSION(:), ALLOCATABLE :: Acoef_U, Bcoef_U + !$OMP THREADPRIVATE(Acoef_U, Bcoef_U) + REAL, DIMENSION(:), ALLOCATABLE :: Acoef_V, Bcoef_V + !$OMP THREADPRIVATE(Acoef_V, Bcoef_V) + LOGICAL :: firstcall=.TRUE. + !$OMP THREADPRIVATE(firstcall) + + + PUBLIC :: climb_wind_down, climb_wind_up + +CONTAINS +! +!**************************************************************************************** +! + SUBROUTINE climb_wind_init + + INTEGER :: ierr + CHARACTER(len = 20) :: modname = 'climb_wind_init' + +!**************************************************************************************** +! Allocation of global module variables +! +!**************************************************************************************** + + ALLOCATE(alf1(klon), stat=ierr) + IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate alf1',1) + + ALLOCATE(alf2(klon), stat=ierr) + IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate alf2',1) + + ALLOCATE(Kcoefm(klon,klev), stat=ierr) + IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate Kcoefm',1) + + ALLOCATE(Ccoef_U(klon,klev), stat=ierr) + IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate Ccoef_U',1) + + ALLOCATE(Dcoef_U(klon,klev), stat=ierr) + IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Dcoef_U',1) + + ALLOCATE(Ccoef_V(klon,klev), stat=ierr) + IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Ccoef_V',1) + + ALLOCATE(Dcoef_V(klon,klev), stat=ierr) + IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Dcoef_V',1) + + ALLOCATE(Acoef_U(klon), Bcoef_U(klon), Acoef_V(klon), Bcoef_V(klon), STAT=ierr) + IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_U and Bcoef_U, ierr=', ierr + + firstcall=.FALSE. + + END SUBROUTINE climb_wind_init +! +!**************************************************************************************** +! + SUBROUTINE climb_wind_down(knon, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old, & +!!! nrlmd le 02/05/2011 + Ccoef_U_out, Ccoef_V_out, Dcoef_U_out, Dcoef_V_out, & + Kcoef_m_out, alf_1_out, alf_2_out, & +!!! + Acoef_U_out, Acoef_V_out, Bcoef_U_out, Bcoef_V_out) +! +! This routine calculates for the wind components u and v, +! recursivly the coefficients C and D in equation +! X(k) = C(k) + D(k)*X(k-1), X=[u,v], k=[1,klev] is the vertical layer. +! +! +USE yomcst_mod_h +USE compbl_mod_h +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: knon + REAL, INTENT(IN) :: dtime + REAL, DIMENSION(klon,klev), INTENT(IN) :: coef_in + REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pres au milieu de couche (Pa) + REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche (Pa) + REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature + REAL, DIMENSION(klon,klev), INTENT(IN) :: delp + REAL, DIMENSION(klon,klev), INTENT(IN) :: u_old + REAL, DIMENSION(klon,klev), INTENT(IN) :: v_old + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_U_out + REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_V_out + REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_U_out + REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_V_out + +!!! nrlmd le 02/05/2011 + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_U_out + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_V_out + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_U_out + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_V_out + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Kcoef_m_out + REAL, DIMENSION(klon), INTENT(OUT) :: alf_1_out + REAL, DIMENSION(klon), INTENT(OUT) :: alf_2_out +!!! + +! Local variables +!**************************************************************************************** + REAL, DIMENSION(klon) :: u1lay, v1lay + INTEGER :: k, i +!**************************************************************************************** +! Initialize module + IF (firstcall) CALL climb_wind_init + +!**************************************************************************************** +! Calculate the coefficients C and D in : u(k) = C(k) + D(k)*u(k-1) +! +!**************************************************************************************** +! - Define alpha (alf1 and alf2) + alf1(:) = 1.0 + alf2(:) = 1.0 - alf1(:) + +! - Calculate the coefficients K + Kcoefm(:,:) = 0.0 + DO k = 2, klev + DO i=1,knon + Kcoefm(i,k) = coef_in(i,k)*RG*RG*dtime/(pplay(i,k-1)-pplay(i,k)) & + *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2 + END DO + END DO + +! - Calculate the coefficients C and D, component "u" + CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), & + u_old(:,:), alf1(:), alf2(:), & + Ccoef_U(:,:), Dcoef_U(:,:), Acoef_U(:), Bcoef_U(:)) + +! - Calculate the coefficients C and D, component "v" + CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), & + v_old(:,:), alf1(:), alf2(:), & + Ccoef_V(:,:), Dcoef_V(:,:), Acoef_V(:), Bcoef_V(:)) + +!**************************************************************************************** +! 6) +! Return the first layer in output variables +! +!**************************************************************************************** + Acoef_U_out = Acoef_U + Bcoef_U_out = Bcoef_U + Acoef_V_out = Acoef_V + Bcoef_V_out = Bcoef_V + +!**************************************************************************************** +! 7) +! If Pbl is split, return also the other layers in output variables +! +!**************************************************************************************** +!!! jyg le 07/02/2012 +!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN + IF (mod(iflag_pbl_split,10) .ge.1) THEN +!!! nrlmd le 02/05/2011 + DO k= 1, klev + DO i= 1, klon + Ccoef_U_out(i,k) = Ccoef_U(i,k) + Ccoef_V_out(i,k) = Ccoef_V(i,k) + Dcoef_U_out(i,k) = Dcoef_U(i,k) + Dcoef_V_out(i,k) = Dcoef_V(i,k) + Kcoef_m_out(i,k) = Kcoefm(i,k) + ENDDO + ENDDO + DO i= 1, klon + alf_1_out(i) = alf1(i) + alf_2_out(i) = alf2(i) + ENDDO +!!! + ENDIF ! (mod(iflag_pbl_split,2) .ge.1) +!!! + + END SUBROUTINE climb_wind_down +! +!**************************************************************************************** +! + SUBROUTINE calc_coef(knon, Kcoef, delp, X, alfa1, alfa2, Ccoef, Dcoef, Acoef, Bcoef) +! +! Find the coefficients C and D in fonction of alfa, K and delp +USE yomcst_mod_h +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: knon + REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp + REAL, DIMENSION(klon,klev), INTENT(IN) :: X + REAL, DIMENSION(klon), INTENT(IN) :: alfa1, alfa2 + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef + +! local variables +!**************************************************************************************** + INTEGER :: k, i + REAL :: buf + !**************************************************************************************** +! + +! Calculate coefficients C and D at top level, k=klev +! + Ccoef(:,:) = 0.0 + Dcoef(:,:) = 0.0 + + DO i = 1, knon + buf = delp(i,klev) + Kcoef(i,klev) + + Ccoef(i,klev) = X(i,klev)*delp(i,klev)/buf + Dcoef(i,klev) = Kcoef(i,klev)/buf + END DO + +! +! Calculate coefficients C and D at top level (klev-1) <= k <= 2 +! + DO k=(klev-1),2,-1 + DO i = 1, knon + buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1)) + + Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1))/buf + Dcoef(i,k) = Kcoef(i,k)/buf + END DO + END DO + +! +! Calculate coeffiecent A and B at surface +! + DO i = 1, knon + buf = delp(i,1) + Kcoef(i,2)*(1-Dcoef(i,2)) + Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*Ccoef(i,2))/buf + Bcoef(i) = -RG/buf + END DO + + END SUBROUTINE calc_coef +! +!**************************************************************************************** +! + + SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, flx_u1, flx_v1, & +!!! nrlmd le 02/05/2011 + Acoef_U_in, Acoef_V_in, Bcoef_U_in, Bcoef_V_in, & + Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in, & + Kcoef_m_in, & +!!! + flx_u_new, flx_v_new, d_u_new, d_v_new) +! +! Diffuse the wind components from the surface layer and up to the top layer. +! Coefficents A, B, C and D are known from before. Start values for the diffusion are the +! momentum fluxes at surface. +! +! u(k=1) = A + B*flx*dtime +! u(k) = C(k) + D(k)*u(k-1) [2 <= k <= klev] +! +!**************************************************************************************** +USE yomcst_mod_h +USE compbl_mod_h +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: knon + REAL, INTENT(IN) :: dtime + REAL, DIMENSION(klon,klev), INTENT(IN) :: u_old + REAL, DIMENSION(klon,klev), INTENT(IN) :: v_old + REAL, DIMENSION(klon), INTENT(IN) :: flx_u1, flx_v1 ! momentum flux + +!!! nrlmd le 02/05/2011 + REAL, DIMENSION(klon), INTENT(IN) :: Acoef_U_in,Acoef_V_in, Bcoef_U_in, Bcoef_V_in + REAL, DIMENSION(klon,klev), INTENT(IN) :: Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in + REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef_m_in +!!! + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon,klev), INTENT(OUT) :: flx_u_new, flx_v_new + REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_u_new, d_v_new + +! Local variables +!**************************************************************************************** + REAL, DIMENSION(klon,klev) :: u_new, v_new + INTEGER :: k, i +!**************************************************************************************** + +!!! jyg le 07/02/2012 +!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN + IF (mod(iflag_pbl_split,10) .ge.1) THEN +!!! nrlmd le 02/05/2011 + DO i = 1, knon + Acoef_U(i)=Acoef_U_in(i) + Acoef_V(i)=Acoef_V_in(i) + Bcoef_U(i)=Bcoef_U_in(i) + Bcoef_V(i)=Bcoef_V_in(i) + ENDDO + DO k = 1, klev + DO i = 1, knon + Ccoef_U(i,k)=Ccoef_U_in(i,k) + Ccoef_V(i,k)=Ccoef_V_in(i,k) + Dcoef_U(i,k)=Dcoef_U_in(i,k) + Dcoef_V(i,k)=Dcoef_V_in(i,k) + Kcoefm(i,k)=Kcoef_m_in(i,k) + ENDDO + ENDDO +!!! + ENDIF ! (mod(iflag_pbl_split,2) .ge.1) +!!! + +! Niveau 1 + DO i = 1, knon + u_new(i,1) = Acoef_U(i) + Bcoef_U(i)*flx_u1(i)*dtime + v_new(i,1) = Acoef_V(i) + Bcoef_V(i)*flx_v1(i)*dtime + END DO + +! Niveau 2 jusqu'au sommet klev + DO k = 2, klev + DO i=1, knon + u_new(i,k) = Ccoef_U(i,k) + Dcoef_U(i,k) * u_new(i,k-1) + v_new(i,k) = Ccoef_V(i,k) + Dcoef_V(i,k) * v_new(i,k-1) + END DO + END DO + +!**************************************************************************************** +! Calcul flux +! +!== flux_u/v est le flux de moment angulaire (positif vers bas) +!== dont l'unite est: (kg m/s)/(m**2 s) +! +!**************************************************************************************** +! + flx_u_new(:,:) = 0.0 + flx_v_new(:,:) = 0.0 + + flx_u_new(1:knon,1)=flx_u1(1:knon) + flx_v_new(1:knon,1)=flx_v1(1:knon) + +! Niveau 2->klev + DO k = 2, klev + DO i = 1, knon + flx_u_new(i,k) = Kcoefm(i,k)/RG/dtime * & + (u_new(i,k)-u_new(i,k-1)) + + flx_v_new(i,k) = Kcoefm(i,k)/RG/dtime * & + (v_new(i,k)-v_new(i,k-1)) + END DO + END DO + +!**************************************************************************************** +! Calcul tendances +! +!**************************************************************************************** + d_u_new(:,:) = 0.0 + d_v_new(:,:) = 0.0 + DO k = 1, klev + DO i = 1, knon + d_u_new(i,k) = u_new(i,k) - u_old(i,k) + d_v_new(i,k) = v_new(i,k) - v_old(i,k) + END DO + END DO + + END SUBROUTINE climb_wind_up +! +!**************************************************************************************** +! +END MODULE climb_wind_mod Index: LMDZ6/trunk/libf/phylmdiso/coef_diff_turb_mod.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/coef_diff_turb_mod.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/coef_diff_turb_mod.f90 (revision 5986) @@ -0,0 +1,603 @@ +! +MODULE coef_diff_turb_mod +! +! This module contains some procedures for calculation of the coefficients of the +! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion +! at surface(cdrag) +! + USE clesphys_mod_h + IMPLICIT NONE + +CONTAINS +! +!**************************************************************************************** +! + SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, & + ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, & + ycoefm, ycoefh ,yq2, yeps, ydrgpro) + + USE dimphy + USE indice_sol_mod + USE print_control_mod, ONLY: prt_level, lunout + USE yomcst_mod_h + USE yoethf_mod_h + USE compbl_mod_h +! +! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the +! atmosphere +! NB! No values are calculated between surface and the first model layer. +! ycoefm(:,1) and ycoefh(:,1) are not valid !!! +! +! +! Input arguments +!**************************************************************************************** + REAL, INTENT(IN) :: dtime + INTEGER, INTENT(IN) :: nsrf, knon + INTEGER, DIMENSION(klon), INTENT(IN) :: ni + REAL, DIMENSION(klon,klev+1), INTENT(IN) :: ypaprs + REAL, DIMENSION(klon,klev), INTENT(IN) :: ypplay + REAL, DIMENSION(klon,klev), INTENT(IN) :: yu, yv + REAL, DIMENSION(klon,klev), INTENT(IN) :: yq, yt + REAL, DIMENSION(klon), INTENT(IN) :: yts, yqsurf + REAL, DIMENSION(klon), INTENT(IN) :: ycdragm +!FC + REAL, DIMENSION(klon,klev), INTENT(IN) :: ydrgpro + + +! InOutput arguments +!**************************************************************************************** + REAL, DIMENSION(klon,klev+1), INTENT(INOUT):: yq2 + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: yeps + REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefh + REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefm + +! Other local variables +!**************************************************************************************** + INTEGER :: k, i, j + REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta + REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq + REAL, DIMENSION(klon) :: yustar + + ykmm = 0 !ym missing init + ykmn = 0 !ym missing init + ykmq = 0 !ym missing init + + yeps(:,:) = 0. + +!**************************************************************************************** +! Calcul de coefficients de diffusion turbulent de l'atmosphere : +! ycoefm(:,2:klev), ycoefh(:,2:klev) +! +!**************************************************************************************** + + CALL coefkz(nsrf, knon, ypaprs, ypplay, & + ksta, ksta_ter, & + yts, yu, yv, yt, yq, & + yqsurf, & + ycoefm, ycoefh) + +!**************************************************************************************** +! Eventuelle recalcule des coeffeicients de diffusion turbulent de l'atmosphere : +! ycoefm(:,2:klev), ycoefh(:,2:klev) +! +!**************************************************************************************** + + IF (iflag_pbl.EQ.1) THEN + CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, & + ycoefm0, ycoefh0) + + DO k = 2, klev + DO i = 1, knon + ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k)) + ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k)) + ENDDO + ENDDO + ENDIF + + +!**************************************************************************************** +! Calcul d'une diffusion minimale pour les conditions tres stables +! +!**************************************************************************************** + IF (ok_kzmin) THEN + CALL coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycdragm, & + ycoefm0,ycoefh0) + + DO k = 2, klev + DO i = 1, knon + ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k)) + ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k)) + ENDDO + ENDDO + + ENDIF + + +!**************************************************************************************** +! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin +! +!**************************************************************************************** + + IF (iflag_pbl.GE.3) THEN + + yzlay(1:knon,1)= & + RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1))) & + *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG + DO k=2,klev + DO i = 1, knon + yzlay(i,k)= & + yzlay(i,k-1)+RD*0.5*(yt(i,k-1)+yt(i,k)) & + /ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG + END DO + END DO + + DO k=1,klev + DO i = 1, knon + yteta(i,k)= & + yt(i,k)*(ypaprs(i,1)/ypplay(i,k))**RKAPPA & + *(1.+0.61*yq(i,k)) + END DO + END DO + + yzlev(1:knon,1)=0. + yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1) + DO k=2,klev + DO i = 1, knon + yzlev(i,k)=0.5*(yzlay(i,k)+yzlay(i,k-1)) + END DO + END DO + +!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!$! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un +!!$! bug sur les coefficients de surface : +!!$! ycdragh(1:knon) = ycoefm(1:knon,1) +!!$! ycdragm(1:knon) = ycoefh(1:knon,1) +!!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! Normalement, on peut passer dans les codes avec knon=0 + ! Mais ca fait planter le replay. + ! En attendant une réécriture, on a joute des if (Fredho) + if ( klon>1 .or. (klon==1 .and. knon==1) ) then + CALL ustarhb(knon,klev,knon,yu,yv,ycdragm, yustar) + endif + + IF (prt_level > 9) THEN + WRITE(lunout,*) 'USTAR = ',(yustar(i),i=1,knon) + ENDIF + +! iflag_pbl peut etre utilise comme longuer de melange + IF (iflag_pbl.GE.31) THEN + if ( klon>1 .or. (klon==1 .and. knon==1) ) then + CALL vdif_kcay(knon,klev,knon,dtime,RG,RD,ypaprs,yt, & + yzlev,yzlay,yu,yv,yteta, & + ycdragm,yq2,q2diag,ykmm,ykmn,yustar, & + iflag_pbl) + endif + ELSE IF (iflag_pbl<20) THEN + CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, & + yzlev,yzlay,yu,yv,yteta, & + ycdragm,yq2,yeps,ykmm,ykmn,ykmq,yustar, & + iflag_pbl,ydrgpro) +!FC + ENDIF + + ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev) + ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev) + + ELSE + ! No TKE for Standard Physics + yq2=0. + ENDIF !(iflag_pbl.ge.3) + + END SUBROUTINE coef_diff_turb +! +!**************************************************************************************** +! + SUBROUTINE coefkz(nsrf, knon, paprs, pplay, & + ksta, ksta_ter, & + ts, & + u,v,t,q, & + qsurf, & + pcfm, pcfh) + + USE yomcst_mod_h + USE dimphy + USE indice_sol_mod + USE print_control_mod, ONLY: prt_level, lunout + USE yoethf_mod_h + USE compbl_mod_h + +!====================================================================== +! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922 +! (une version strictement identique a l'ancien modele) +! Objet: calculer le coefficient du frottement du sol (Cdrag) et les +! coefficients d'echange turbulent dans l'atmosphere. +! Arguments: +! nsrf-----input-I- indicateur de la nature du sol +! knon-----input-I- nombre de points a traiter +! paprs----input-R- pregssion a chaque intercouche (en Pa) +! pplay----input-R- pression au milieu de chaque couche (en Pa) +! ts-------input-R- temperature du sol (en Kelvin) +! u--------input-R- vitesse u +! v--------input-R- vitesse v +! t--------input-R- temperature (K) +! q--------input-R- vapeur d'eau (kg/kg) +! +! pcfm-----output-R- coefficients a calculer (vitesse) +! pcfh-----output-R- coefficients a calculer (chaleur et humidite) +!====================================================================== + INCLUDE "FCTTRE.h" +! +! Arguments: +! + INTEGER, INTENT(IN) :: knon, nsrf + REAL, INTENT(IN) :: ksta, ksta_ter + REAL, DIMENSION(klon), INTENT(IN) :: ts + REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs + REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay + REAL, DIMENSION(klon,klev), INTENT(IN) :: u, v, t, q + REAL, DIMENSION(klon), INTENT(IN) :: qsurf + + REAL, DIMENSION(klon,klev), INTENT(OUT) :: pcfm, pcfh + +! +! Local variables: +! + INTEGER, DIMENSION(klon) :: itop ! numero de couche du sommet de la couche limite +! +! Quelques constantes et options: +! + REAL, PARAMETER :: cepdu2=0.1**2 + REAL, PARAMETER :: CKAP=0.4 + REAL, PARAMETER :: cb=5.0 + REAL, PARAMETER :: cc=5.0 + REAL, PARAMETER :: cd=5.0 + REAL, PARAMETER :: clam=160.0 + REAL, PARAMETER :: ratqs=0.05 ! largeur de distribution de vapeur d'eau + LOGICAL, PARAMETER :: richum=.TRUE. ! utilise le nombre de Richardson humide + REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique + REAL, PARAMETER :: prandtl=0.4 + REAL kstable ! diffusion minimale (situation stable) + ! GKtest + ! PARAMETER (kstable=1.0e-10) +!IM: 261103 REAL kstable_ter, kstable_sinon +!IM: 211003 cf GK PARAMETER (kstable_ter = 1.0e-6) +!IM: 261103 PARAMETER (kstable_ter = 1.0e-8) +!IM: 261103 PARAMETER (kstable_ter = 1.0e-10) +!IM: 261103 PARAMETER (kstable_sinon = 1.0e-10) + ! fin GKtest + REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange + INTEGER isommet ! le sommet de la couche limite + LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante + LOGICAL, PARAMETER :: opt_ec=.FALSE.! formule du Centre Europeen dans l'atmosphere + +! +! Variables locales: + INTEGER i, k !IM 120704 + REAL zgeop(klon,klev) + REAL zmgeom(klon) + REAL zri(klon) + REAL zl2(klon) + REAL zdphi, zdu2, ztvd, ztvu, zcdn + REAL zscf + REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs + REAL z2geomf, zalh2, zalm2, zscfh, zscfm + REAL, PARAMETER :: t_coup=273.15 + LOGICAL, PARAMETER :: check=.FALSE. +! +! contre-gradient pour la chaleur sensible: Kelvin/metre + REAL gamt(2:klev) + + LOGICAL, SAVE :: appel1er=.TRUE. + !$OMP THREADPRIVATE(appel1er) +! +! Fonctions thermodynamiques et fonctions d'instabilite + REAL fsta, fins, x + + fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x)) + fins(x) = SQRT(1.0-18.0*x) + + isommet=klev + + IF (appel1er) THEN + IF (prt_level > 9) THEN + WRITE(lunout,*)'coefkz, opt_ec:', opt_ec + WRITE(lunout,*)'coefkz, richum:', richum + IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs + WRITE(lunout,*)'coefkz, isommet:', isommet + WRITE(lunout,*)'coefkz, tvirtu:', tvirtu + appel1er = .FALSE. + ENDIF + ENDIF +! +! Initialiser les sorties +! + DO k = 1, klev + DO i = 1, knon + pcfm(i,k) = 0.0 + pcfh(i,k) = 0.0 + ENDDO + ENDDO + DO i = 1, knon + itop(i) = 0 + ENDDO + +! +! Prescrire la valeur de contre-gradient +! + IF (iflag_pbl.EQ.1) THEN + DO k = 3, klev + gamt(k) = -1.0E-03 + ENDDO + gamt(2) = -2.5E-03 + ELSE + DO k = 2, klev + gamt(k) = 0.0 + ENDDO + ENDIF +!IM cf JLD/ GKtest + IF ( nsrf .NE. is_oce ) THEN +!IM 261103 kstable = kstable_ter + kstable = ksta_ter + ELSE +!IM 261103 kstable = kstable_sinon + kstable = ksta + ENDIF +!IM cf JLD/ GKtest fin + +! +! Calculer les geopotentiels de chaque couche +! + DO i = 1, knon + zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) & + * (paprs(i,1)-pplay(i,1)) + ENDDO + DO k = 2, klev + DO i = 1, knon + zgeop(i,k) = zgeop(i,k-1) & + + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) & + * (pplay(i,k-1)-pplay(i,k)) + ENDDO + ENDDO + +! +! Calculer les coefficients turbulents dans l'atmosphere +! + DO i = 1, knon + itop(i) = isommet + ENDDO + + + DO k = 2, isommet + DO i = 1, knon + zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 & + +(v(i,k)-v(i,k-1))**2) + zmgeom(i)=zgeop(i,k)-zgeop(i,k-1) + zdphi =zmgeom(i) / 2.0 + zt = (t(i,k)+t(i,k-1)) * 0.5 + zq = (q(i,k)+q(i,k-1)) * 0.5 + +! +! Calculer Qs et dQs/dT: +! + IF (thermcep) THEN + zdelta = MAX(0.,SIGN(1.,RTT-zt)) + zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) & + + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta + zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k) + zqs = MIN(0.5,zqs) + zcor = 1./(1.-RETV*zqs) + zqs = zqs*zcor + zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor) + ELSE + IF (zt .LT. t_coup) THEN + zqs = qsats(zt) / pplay(i,k) + zdqs = dqsats(zt,zqs) + ELSE + zqs = qsatl(zt) / pplay(i,k) + zdqs = dqsatl(zt,zqs) + ENDIF + ENDIF +! +! calculer la fraction nuageuse (processus humide): +! + if (zq /= 0.) then + zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq) + else + zfr = 0. + end if + zfr = MAX(0.0,MIN(1.0,zfr)) + IF (.NOT.richum) zfr = 0.0 +! +! calculer le nombre de Richardson: +! + IF (tvirtu) THEN + ztvd =( t(i,k) & + + zdphi/RCPD/(1.+RVTMP2*zq) & + *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) & + )*(1.+RETV*q(i,k)) + ztvu =( t(i,k-1) & + - zdphi/RCPD/(1.+RVTMP2*zq) & + *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) & + )*(1.+RETV*q(i,k-1)) + zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu)) + zri(i) = zri(i) & + + zmgeom(i)*zmgeom(i)/RG*gamt(k) & + *(paprs(i,k)/101325.0)**RKAPPA & + /(zdu2*0.5*(ztvd+ztvu)) + + ELSE ! calcul de Ridchardson compatible LMD5 + + zri(i) =(RCPD*(t(i,k)-t(i,k-1)) & + -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k) & + *(pplay(i,k)-pplay(i,k-1)) & + )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k))) + zri(i) = zri(i) + & + zmgeom(i)*zmgeom(i)*gamt(k)/RG & + *(paprs(i,k)/101325.0)**RKAPPA & + /(zdu2*0.5*(t(i,k-1)+t(i,k))) + ENDIF +! +! finalement, les coefficients d'echange sont obtenus: +! + zcdn=SQRT(zdu2) / zmgeom(i) * RG + + IF (opt_ec) THEN + z2geomf=zgeop(i,k-1)+zgeop(i,k) + zalm2=(0.5*ckap/RG*z2geomf & + /(1.+0.5*ckap/rg/clam*z2geomf))**2 + zalh2=(0.5*ckap/rg*z2geomf & + /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2 + IF (zri(i).LT.0.0) THEN ! situation instable + zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3 & + / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG) + zscf = SQRT(-zri(i)*zscf) + zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf) + zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf) + pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm) + pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh) + ELSE ! situation stable + zscf=SQRT(1.+cd*zri(i)) + pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf) + pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf) + ENDIF + ELSE + zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) & + /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2 + pcfm(i,k)=SQRT(MAX(zcdn*zcdn*(ric-zri(i))/ric, kstable)) + pcfm(i,k)= zl2(i)* pcfm(i,k) + pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different + ENDIF + ENDDO + ENDDO + +! +! Au-dela du sommet, pas de diffusion turbulente: +! + DO i = 1, knon + IF (itop(i)+1 .LE. klev) THEN + DO k = itop(i)+1, klev + pcfh(i,k) = 0.0 + pcfm(i,k) = 0.0 + ENDDO + ENDIF + ENDDO + + END SUBROUTINE coefkz +! +!**************************************************************************************** +! + SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, & + pcfm, pcfh) + + USE yomcst_mod_h +USE dimphy + USE indice_sol_mod + +!====================================================================== +! J'introduit un peu de diffusion sauf dans les endroits +! ou une forte inversion est presente +! On peut dire qu'il represente la convection peu profonde +! +! Arguments: +! nsrf-----input-I- indicateur de la nature du sol +! knon-----input-I- nombre de points a traiter +! paprs----input-R- pression a chaque intercouche (en Pa) +! pplay----input-R- pression au milieu de chaque couche (en Pa) +! t--------input-R- temperature (K) +! +! pcfm-----output-R- coefficients a calculer (vitesse) +! pcfh-----output-R- coefficients a calculer (chaleur et humidite) +!====================================================================== +! +! Arguments: +! + INTEGER, INTENT(IN) :: knon, nsrf + REAL, DIMENSION(klon, klev+1), INTENT(IN) :: paprs + REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay + REAL, DIMENSION(klon, klev), INTENT(IN) :: t(klon,klev) + + REAL, DIMENSION(klon, klev), INTENT(OUT) :: pcfm, pcfh +! +! Quelques constantes et options: +! + REAL, PARAMETER :: prandtl=0.4 + REAL, PARAMETER :: kstable=0.002 +! REAL, PARAMETER :: kstable=0.001 + REAL, PARAMETER :: mixlen=35.0 ! constante controlant longueur de melange + REAL, PARAMETER :: seuil=-0.02 ! au-dela l'inversion est consideree trop faible +! PARAMETER (seuil=-0.04) +! PARAMETER (seuil=-0.06) +! PARAMETER (seuil=-0.09) + +! +! Variables locales: +! + INTEGER i, k, invb(knon) + REAL zl2(knon) + REAL zdthmin(knon), zdthdp + + +! +! Initialiser les sorties +! + DO k = 1, klev + DO i = 1, knon + pcfm(i,k) = 0.0 + pcfh(i,k) = 0.0 + ENDDO + ENDDO + +! +! Chercher la zone d'inversion forte +! + DO i = 1, knon + invb(i) = klev + zdthmin(i)=0.0 + ENDDO + DO k = 2, klev/2-1 + DO i = 1, knon + zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) & + - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1) + zdthdp = zdthdp * 100.0 + IF (pplay(i,k).GT.0.8*paprs(i,1) .AND. & + zdthdp.LT.zdthmin(i) ) THEN + zdthmin(i) = zdthdp + invb(i) = k + ENDIF + ENDDO + ENDDO + +! +! Introduire une diffusion: +! + IF ( nsrf.EQ.is_oce ) THEN + DO k = 2, klev + DO i = 1, knon +!IM cf FH/GK IF ( (nsrf.NE.is_oce) .OR. ! si ce n'est pas sur l'ocean +!IM cf FH/GK . (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion + !IM cf JLD/ GKtest TERkz2 + ! IF ( (nsrf.EQ.is_ter) .OR. ! si on est sur la terre + ! fin GKtest + + +! s'il n'y a pas d'inversion ou si l'inversion est trop faible +! IF ( (nsrf.EQ.is_oce) .AND. & + IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN + zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1)) & + /(paprs(i,2)-paprs(i,klev+1)) ))**2 + pcfm(i,k)= zl2(i)* kstable + pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different + ENDIF + ENDDO + ENDDO + ENDIF + + END SUBROUTINE coefkz2 +! +!**************************************************************************************** +! +END MODULE coef_diff_turb_mod Index: LMDZ6/trunk/libf/phylmdiso/coefkzmin.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/coefkzmin.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/coefkzmin.f90 (revision 5986) @@ -0,0 +1,133 @@ + +SUBROUTINE coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycdragm, km, kn) + + USE dimphy + USE yomcst_mod_h +IMPLICIT NONE + + + + ! ....................................................................... + ! Entrees modifies en attendant une version ou les zlev, et zlay soient + ! disponibles. + + REAL ycdragm(klon) + + REAL yu(klon, klev), yv(klon, klev) + REAL yt(klon, klev), yq(klon, klev) + REAL ypaprs(klon, klev+1), ypplay(klon, klev) + REAL yustar(klon) + REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev) + + INTEGER i + + ! ....................................................................... + + ! En entree : + ! ----------- + + ! zlev : altitude a chaque niveau (interface inferieure de la couche + ! de meme indice) + ! ustar : u* + + ! teta : temperature potentielle au centre de chaque couche + ! (en entree : la valeur au debut du pas de temps) + + ! en sortier : + ! ------------ + + ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque + ! couche) + ! (en sortie : la valeur a la fin du pas de temps) + ! kn : diffusivite turbulente des scalaires (au bas de chaque couche) + ! (en sortie : la valeur a la fin du pas de temps) + + ! ....................................................................... + + REAL ustar(klon) + REAL kmin, qmin, pblhmin(klon), coriol(klon) + REAL zlev(klon, klev+1) + REAL teta(klon, klev) + + REAL km(klon, klev) + REAL kn(klon, klev) + INTEGER knon + + + INTEGER nlay, nlev + INTEGER ig, k + + REAL, PARAMETER :: kap = 0.4 + + nlay = klev + nlev = klev + 1 + ! ....................................................................... + ! en attendant une version ou les zlev, et zlay soient + ! disponibles. + ! Debut de la partie qui doit etre unclue a terme dans clmain. + + DO i = 1, knon + yzlay(i, 1) = rd*yt(i, 1)/(0.5*(ypaprs(i,1)+ypplay(i, & + 1)))*(ypaprs(i,1)-ypplay(i,1))/rg + END DO + DO k = 2, klev + DO i = 1, knon + yzlay(i, k) = yzlay(i, k-1) + rd*0.5*(yt(i,k-1)+yt(i,k))/ypaprs(i, k)*( & + ypplay(i,k-1)-ypplay(i,k))/rg + END DO + END DO + DO k = 1, klev + DO i = 1, knon + ! ATTENTION:on passe la temperature potentielle virt. pour le calcul de + ! K + yteta(i, k) = yt(i, k)*(ypaprs(i,1)/ypplay(i,k))**rkappa* & + (1.+0.61*yq(i,k)) + END DO + END DO + DO i = 1, knon + yzlev(i, 1) = 0. + yzlev(i, klev+1) = 2.*yzlay(i, klev) - yzlay(i, klev-1) + END DO + DO k = 2, klev + DO i = 1, knon + yzlev(i, k) = 0.5*(yzlay(i,k)+yzlay(i,k-1)) + END DO + END DO + + yustar(1:knon) = sqrt(ycdragm(1:knon)*(yu(1:knon,1)*yu(1:knon,1)+yv(1:knon, & + 1)*yv(1:knon,1))) + + ! Fin de la partie qui doit etre unclue a terme dans clmain. + + ! ette routine est ecrite pour avoir en entree ustar, teta et zlev + ! Ici, on a inclut le calcul de ces trois variables dans la routine + ! coefkzmin en attendant une nouvelle version de la couche limite + ! ou ces variables seront disponibles. + + ! Debut de la routine coefkzmin proprement dite. + + ustar = yustar + teta = yteta + zlev = yzlev + + DO ig = 1, knon + coriol(ig) = 1.E-4 + pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5) + END DO + + DO k = 2, klev + DO ig = 1, knon + IF (teta(ig,2)>teta(ig,1)) THEN + qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2 + kmin = kap*zlev(ig, k)*qmin + ELSE + kmin = 0. ! kmin n'est utilise que pour les SL stables. + END IF + kn(ig, k) = kmin + km(ig, k) = kmin + END DO + END DO + + + RETURN +END SUBROUTINE coefkzmin Index: LMDZ6/trunk/libf/phylmdiso/freinage.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/freinage.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/freinage.f90 (revision 5986) @@ -0,0 +1,135 @@ +! +! $Id$ +! + SUBROUTINE freinage(knon, uu, vv, & + tt,veget,lai, height,ypaprs,ypplay,drag_pro,d_u,d_v) + + !ONLINE: +USE yoegwd_mod_h + USE dimpft_mod_h + USE compbl_mod_h + USE clesphys_mod_h + use dimphy, only: klon, klev +! USE control, ONLY: nvm +! USE indice_sol_mod, only : nvm_orch + + USE yomcst_mod_h +IMPLICIT NONE + + + +!FC + + ! 0. DECLARATIONS: + + ! 0.1 INPUTS + + REAL, DIMENSION(klon,klev), INTENT(IN) :: ypplay + REAL, DIMENSION(klon,klev+1), INTENT(IN) :: ypaprs + + + REAL, DIMENSION(klon, klev), INTENT(IN) :: uu + REAL, DIMENSION(klon, klev), INTENT(IN) :: vv + REAL, DIMENSION(klon, klev), INTENT(IN) :: tt + REAL, DIMENSION(klon,nvm_lmdz), INTENT(IN) :: veget,lai + REAL, DIMENSION(klon,nvm_lmdz), INTENT(IN) :: height + + REAL, DIMENSION(klon,klev) :: wind + REAL, DIMENSION(klon, klev) :: yzlay + INTEGER knon + + ! 0.2 OUTPUTS + + REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v + REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in v + !knon nombre de points concernes + REAL, DIMENSION(klon,klev) :: sumveg ! change in v + + REAL, DIMENSION(klon,klev), INTENT(OUT) :: drag_pro + ! (KLON, KLEV) tendencies on winds + + + INTEGER k,jv,i + + +!FCCCC REAL Cd_frein + + ! 0.3.1 LOCAL VARIABLE + + + !----------------------------------------------------------------- + + ! 1. INITIALISATIONS + + +! Cd_frein = 7.5E-2 ! (0.075) ! Drag from MASSON 2009 +!FC ESSAI +! Cd_frein = 1.5E-2 ! (0.075) ! Drag from MASSON 2009 +! Cd_frein = 0.005 ! (0.075) ! Drag from MASSON 2009 + +! initialisation + d_u(:,:) =0. + d_v(:,:) =0. + drag_pro(:,:) =0. + sumveg(:,:) =0. +!! print*, "Cd_frein" , Cd_frein + + wind(:,:)= sqrt(uu(:,:)*uu(:,:)+vv(:,:)*vv(:,:)) + + yzlay(1:knon,1)= & + RD*tt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1))) & + *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG + DO k=2,klev + yzlay(1:knon,k)= & + yzlay(1:knon,k-1)+RD*0.5*(tt(1:knon,k-1)+tt(1:knon,k)) & + /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG + END DO + +! verifier les indexes ..... +!! print*, " calcul de drag_pro FC " + + do k= 1,klev + + do jv=2,nvm_lmdz ! (on peut faire 9 ?) + + do i=1,knon + + sumveg(i,k)= sumveg(i,k)+ veget(i,jv) + +! if ( (height(i,jv) .gt. yzlay(i,k)) .AND. (height(i,jv) .gt. 0.1) .and. LAI(i,jv).gt.0. ) then + if ( (height(i,jv) .gt. yzlay(i,k)) .AND. (height(i,jv) .gt. 0.1) ) then +!FC attention veut on le test sur le LAI ? + if (ifl_pbltree.eq.1) then + drag_pro(i,k)= drag_pro(i,k)+ & + veget(i,jv) + elseif (ifl_pbltree.eq.2) then + drag_pro(i,k)= drag_pro(i,k)+ & + 6*LAI(i,jv)*veget(i,jv)*( yzlay(i,k)*(height(i,jv)-yzlay(i,k))/(height(i,jv)*height(i,jv)+ 0.01)) + elseif (ifl_pbltree.eq.3) then + drag_pro(i,k)= drag_pro(i,k)+ & + veget(i,jv)*( yzlay(i,k)*(height(i,jv)-yzlay(i,k))/(height(i,jv)*height(i,jv)+ 0.01)) + elseif (ifl_pbltree.eq.0) then + drag_pro(i,k)=0.0 + endif + else + drag_pro(i,k)= drag_pro(i,k) + endif + + + enddo + enddo + enddo + do k=1,klev + where (sumveg(1:knon,k) > 0.05 ) +! drag_pro(1:knon,k)=Cd_frein*drag_pro(1:knon,k)/sumveg(1:knon,k) + drag_pro(1:knon,k)=Cd_frein*drag_pro(1:knon,k) + elsewhere + drag_pro(1:knon,k)=0.0 + endwhere + d_u(1:knon,k) =(-1)*drag_pro(1:knon,k)*uu(1:knon,k)*wind(1:knon,k) + d_v(1:knon,k) =(-1)*drag_pro(1:knon,k)*vv(1:knon,k)*wind(1:knon,k) + enddo + return + + END SUBROUTINE freinage + Index: LMDZ6/trunk/libf/phylmdiso/ocean_albedo.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/ocean_albedo.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/ocean_albedo.f90 (revision 5986) @@ -0,0 +1,256 @@ +! +! $Id$ +! + +SUBROUTINE ocean_albedo(knon,zrmu0,knindex,pwind,SFRWL,alb_dir_new,alb_dif_new) +!! +!!**** *ALBEDO_RS14* +!! +!! PURPOSE +!! ------- +!! computes the direct & diffuse albedo over open water +!! +!!** METHOD +!! ------ +!! +!! EXTERNAL +!! -------- +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! R. Séférian * Meteo-France * +!! +!! MODIFICATIONS +!! ------------- +!! Original 03/2014 +!! 05/2014 R. Séférian & B. Decharme :: Adaptation to spectral +!! computation for diffuse and direct albedo +!! 08/2014 S. Baek :: for wider wavelength range 200-4000nm and +!! adaptation to LMDZ + whitecap effect by Koepke + chrolophyll +!! map from climatology file +!! 10/2016 O. Boucher :: some optimisation following R. +!! Seferian's work in the CNRM Model +!! +!------------------------------------------------------------------------------- +! +!* DECLARATIONS +! ------------ +! +USE ocean_albedo_para +USE dimphy +USE phys_state_var_mod, ONLY : chl_con +USE clesphys_mod_h +! +! +IMPLICIT NONE +! +!* 0.1 declarations of arguments +! ------------------------- +! +! +INTEGER, INTENT(IN) :: knon +INTEGER, DIMENSION(klon), INTENT(IN) :: knindex +REAL, DIMENSION(klon), INTENT(IN) :: zrmu0 !--cos(SZA) on full vector +REAL, DIMENSION(klon), INTENT(IN) :: pwind !--wind speed on compressed vector +REAL, DIMENSION(6),INTENT(IN) :: SFRWL +REAL, DIMENSION(klon,nsw), INTENT(OUT) :: alb_dir_new, alb_dif_new +! +!* 0.2 declarations of local variables +! ------------------------- +! +REAL, DIMENSION(klon) :: ZCHL ! surface chlorophyll +REAL, DIMENSION(klon) :: ZCOSZEN ! Cosine of the zenith solar angle +! +INTEGER :: JWL, INU ! indexes +INTEGER :: JI +REAL :: ZWL ! input parameter: wavelength and diffuse/direct fraction of light +REAL:: ZCHLABS, ZAW, ZBW, ZREFM, ZYLMD, ZUE, ZUE2 ! scalar computation variables +! +REAL, DIMENSION(klon) :: ZAP, ZXX2, ZR00, ZRR0, ZRRR ! computation variables +REAL, DIMENSION(klon) :: ZR22, ZR11DF ! computation variables +REAL, DIMENSION(klon) :: ZBBP, ZNU, ZHB ! computation variables +REAL, DIMENSION(klon) :: ZR11, ZRW, ZRWDF, ZRDF ! 4 components of the OSA +REAL, DIMENSION(klon) :: ZSIG, ZFWC, ZWORK1, ZWORK2, ZWORK3 +! +!--initialisations------------- +! + +IF (knon==0) RETURN ! A verifier pourquoi on en a besoin... + +alb_dir_new(:,:) = 0. +alb_dif_new(:,:) = 0. +! +! Initialisation of chlorophyll content +! ZCHL(:) = CHL_CON!0.05 ! averaged global values for surface chlorophyll +IF (ok_chlorophyll) THEN + ZCHL(1:knon)=CHL_CON(knindex(1:knon)) +ELSE + ZCHL(1:knon) = 0.05 +ENDIF + +! variables that do not depend on wavelengths +! loop over the grid points +! functions of chlorophyll content +ZWORK1(1:knon)= EXP(LOG(ZCHL(1:knon))*0.65) +ZWORK2(1:knon)= 0.416 * EXP(LOG(ZCHL(1:knon))*0.766) +ZWORK3(1:knon)= LOG10(ZCHL(1:knon)) +! store the cosine of the solar zenith angle +ZCOSZEN(1:knon) = zrmu0(knindex(1:knon)) +! Compute sigma derived from wind speed (Cox & Munk reflectance model) +ZSIG(1:knon)=SQRT(0.003+0.00512*PWIND(1:knon)) +! original : correction for foam (Eq 16-17) +! has to be update once we have information from wave model (discussion with G. Madec) +ZFWC(1:knon)=3.97e-4*PWIND(1:knon)**1.59 ! Salisbury 2014 eq(2) at 37GHz, value in fraction +! +DO JWL=1,NNWL ! loop over the wavelengths + ! + !--------------------------------------------------------------------------------- + ! 0- Compute baseline values + !--------------------------------------------------------------------------------- + + ! Get refractive index for the correspoding wavelength + ZWL=XAKWL(JWL) !!!--------- wavelength value + ZREFM= XAKREFM(JWL) !!!--------- refraction index value + + !--------------------------------------------------------------------------------- + ! 1- Compute direct surface albedo (ZR11) + !--------------------------------------------------------------------------------- + ! + ZXX2(1:knon)=SQRT(1.0-(1.0-ZCOSZEN(1:knon)**2)/ZREFM**2) + ZRR0(1:knon)=0.50*(((ZXX2(1:knon)-ZREFM*ZCOSZEN(1:knon))/(ZXX2(1:knon)+ZREFM*ZCOSZEN(1:knon)))**2 + & + ((ZCOSZEN(1:knon)-ZREFM*ZXX2(1:knon))/(ZCOSZEN(1:knon)+ZREFM*ZXX2(1:knon)))**2) + ZRRR(1:knon)=0.50*(((ZXX2(1:knon)-1.34*ZCOSZEN(1:knon))/(ZXX2(1:knon)+1.34*ZCOSZEN(1:knon)))**2 + & + ((ZCOSZEN(1:knon)-1.34*ZXX2(1:knon))/(ZCOSZEN(1:knon)+1.34*ZXX2(1:knon)))**2) + ZR11(1:knon)=ZRR0(1:knon)-(0.0152-1.7873*ZCOSZEN(1:knon)+6.8972*ZCOSZEN(1:knon)**2-8.5778*ZCOSZEN(1:knon)**3+ & + 4.071*ZSIG(1:knon)-7.6446*ZCOSZEN(1:knon)*ZSIG(1:knon)) * & + EXP(0.1643-7.8409*ZCOSZEN(1:knon)-3.5639*ZCOSZEN(1:knon)**2-2.3588*ZSIG(1:knon)+ & + 10.0538*ZCOSZEN(1:knon)*ZSIG(1:knon))*ZRR0(1:knon)/ZRRR(1:knon) + ! + !--------------------------------------------------------------------------------- + ! 2- Compute surface diffuse albedo (ZRDF) + !--------------------------------------------------------------------------------- + ! Diffuse albedo from Jin et al., 2006 + estimation from diffuse fraction of + ! light (relying later on AOD). CNRM model has opted for Eq 5b + ZRDF(1:knon)=-0.1482-0.012*ZSIG(1:knon)+0.1609*ZREFM-0.0244*ZSIG(1:knon)*ZREFM ! surface diffuse (Eq 5a) + !!ZRDF(1:knon)=-0.1479+0.1502*ZREFM-0.0176*ZSIG(1:knon)*ZREFM ! surface diffuse (Eq 5b) + + !--------------------------------------------------------------------------------- + ! *- Determine absorption and backscattering + ! coefficients to determine reflectance below the surface (Ro) once for all + ! + ! *.1- Absorption by chlorophyll + ZCHLABS= XAKACHL(JWL) + ! *.2- Absorption by seawater + ZAW= XAKAW3(JWL) + ! *.3- Backscattering by seawater + ZBW= XAKBW(JWL) + ! *.4- Backscattering by chlorophyll + ZYLMD = EXP(0.014*(440.0-ZWL)) + ZAP(1:knon) = 0.06*ZCHLABS*ZWORK1(1:knon) +0.2*(XAW440+0.06*ZWORK1(1:knon))*ZYLMD + +!! WHERE ( ZCHL(1:knon) > 0.02 ) +!! ZNU(:)=MIN(0.0,0.5*(ZWORK3(:)-0.3)) +!! ZBBP(:)=(0.002+0.01*(0.5-0.25*ZWORK3(:))*(ZWL/550.)**ZNU(:))*ZWORK2(:) +!! ELSEWHERE +!! ZBBP(:)=0.019*(550./ZWL)*ZWORK2(:) !ZBBPf=0.0113 at chl<=0.02 +!! ENDWHERE + + do JI = 1, knon + IF (ZCHL(JI) > 0.02) THEN + ZNU(JI)=MIN(0.0,0.5*(ZWORK3(JI)-0.3)) + ZBBP(JI)=(0.002+0.01*(0.5-0.25*ZWORK3(JI))*(ZWL/550.)**ZNU(JI)) & + *ZWORK2(JI) + ELSE + ZBBP(JI)=0.019*(550./ZWL)*ZWORK2(JI) !ZBBPf=0.0113 at chl<=0.02 + ENDIF + ENDDO + + ! Morel-Gentili(1991), Eq (12) + ! ZHB=h/(h+2*ZBBPf*(1.-h)) + ZHB(1:knon)=0.5*ZBW/(0.5*ZBW+ZBBP(1:knon)) + + !--------------------------------------------------------------------------------- + ! 3- Compute direct water-leaving albedo (ZRW) + !--------------------------------------------------------------------------------- + ! Based on Morel & Gentilli 1991 parametrization + ZR22(1:knon)=0.48168549-0.014894708*ZSIG(1:knon)-0.20703885*ZSIG(1:knon)**2 + + ! Use Morel 91 formula to compute the direct reflectance + ! below the surface + ZR00(1:knon)=(0.5*ZBW+ZBBP(1:knon))/(ZAW+ZAP(1:knon)) * & + (0.6279-0.2227*ZHB(1:knon)-0.0513*ZHB(1:knon)**2 + & + (-0.3119+0.2465*ZHB(1:knon))*ZCOSZEN(1:knon)) + ZRW(1:knon)=ZR00(1:knon)*(1.-ZR22(1:knon))/(1.-ZR00(1:knon)*ZR22(1:knon)) + + !--------------------------------------------------------------------------------- + ! 4- Compute diffuse water-leaving albedo (ZRWDF) + !--------------------------------------------------------------------------------- + ! as previous water-leaving computation but assumes a uniform incidence of + ! shortwave at surface (ue) + ZUE=0.676 ! equivalent u_unif for diffuse incidence + ZUE2=SQRT(1.0-(1.0-ZUE**2)/ZREFM**2) + ZRR0(1:knon)=0.50*(((ZUE2-ZREFM*ZUE)/(ZUE2+ZREFM*ZUE))**2 +((ZUE-ZREFM*ZUE2)/(ZUE+ZREFM*ZUE2))**2) + ZRRR(1:knon)=0.50*(((ZUE2-1.34*ZUE)/(ZUE2+1.34*ZUE))**2 +((ZUE-1.34*ZUE2)/(ZUE+1.34*ZUE2))**2) + ZR11DF(1:knon)=ZRR0(1:knon)-(0.0152-1.7873*ZUE+6.8972*ZUE**2-8.5778*ZUE**3+4.071*ZSIG(1:knon)-7.6446*ZUE*ZSIG(1:knon)) * & + EXP(0.1643-7.8409*ZUE-3.5639*ZUE**2-2.3588*ZSIG(1:knon)+10.0538*ZUE*ZSIG(1:knon))*ZRR0(1:knon)/ZRRR(1:knon) + + ! Use Morel 91 formula to compute the diffuse + ! reflectance below the surface + ZR00(1:knon) = (0.5*ZBW+ZBBP(1:knon)) / (ZAW+ZAP(1:knon)) & + * (0.6279-0.2227*ZHB(1:knon)-0.0513*ZHB(1:knon)**2 & + + (-0.3119+0.2465*ZHB(1:knon))*ZUE) + ZRWDF(1:knon)=ZR00(1:knon)*(1.-ZR22(1:knon))*(1.-ZR11DF(1:knon))/(1.-ZR00(1:knon)*ZR22(1:knon)) + + ! get waveband index inu for each nsw band + SELECT CASE(nsw) + CASE(2) + IF (JWL.LE.49) THEN ! from 200 to 680 nm + inu=1 + ELSE ! from 690 to 4000 nm + inu=2 + ENDIF + CASE(4) + IF (JWL.LE.49) THEN ! from 200 to 680 nm + inu=1 + ELSE IF (JWL.LE.99) THEN ! from 690 to 1180 nm + inu=2 + ELSE IF (JWL.LE.218) THEN ! from 1190 to 2370 nm + inu=3 + ELSE ! from 2380 to 4000 nm + inu=4 + ENDIF + CASE(6) + IF (JWL.LE.5) THEN ! from 200 to 240 nm + inu=1 + ELSE IF (JWL.LE.24) THEN ! from 250 to 430 nm + inu=2 + ELSE IF (JWL.LE.49) THEN ! from 440 to 680 nm + inu=3 + ELSE IF (JWL.LE.99) THEN ! from 690 to 1180 nm + inu=4 + ELSE IF (JWL.LE.218) THEN ! from 1190 to 2370 nm + inu=5 + ELSE ! from 2380 to 4000 nm + inu=6 + ENDIF + END SELECT + + ! partitionning direct and diffuse albedo + ! excluding diffuse albedo ZRW on ZDIR_ALB + + !--direct + alb_dir_new(1:knon,inu)=alb_dir_new(1:knon,inu) + & + ( XFRWL(JWL) * ((1.-ZFWC(1:knon)) * (ZR11(1:knon)+ZRW(1:knon)) + ZFWC(1:knon)*XRWC(JWL)) )/SFRWL(inu) + !--diffuse + alb_dif_new(1:knon,inu)=alb_dif_new(1:knon,inu) + & + ( XFRWL(JWL) * ((1.-ZFWC(1:knon)) * (ZRDF(1:knon)+ZRWDF(1:knon)) + ZFWC(1:knon)*XRWC(JWL)) )/SFRWL(inu) + +ENDDO ! ending loop over wavelengths + +END SUBROUTINE ocean_albedo Index: LMDZ6/trunk/libf/phylmdiso/ocean_slab_mod.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/ocean_slab_mod.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/ocean_slab_mod.f90 (revision 5986) @@ -0,0 +1,1018 @@ +!Completed +MODULE ocean_slab_mod +! +! This module is used for both surface ocean and sea-ice when using the slab ocean, +! "ocean=slab". +! + + USE dimphy + USE indice_sol_mod + USE surface_data + USE mod_grid_phy_lmdz, ONLY: klon_glo + USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root + + IMPLICIT NONE + PRIVATE + PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice + +!*********************************************************************************** +! Global saved variables +!*********************************************************************************** + ! number of slab vertical layers + INTEGER, PUBLIC, SAVE :: nslay + !$OMP THREADPRIVATE(nslay) + ! timestep for coupling (update slab temperature) in timesteps + INTEGER, PRIVATE, SAVE :: cpl_pas + !$OMP THREADPRIVATE(cpl_pas) + ! cyang = 1/heat capacity of top layer (rho.c.H) + REAL, PRIVATE, SAVE :: cyang + !$OMP THREADPRIVATE(cyang) + ! depth of slab layers (1 or 2) + REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh + !$OMP THREADPRIVATE(slabh) + ! slab temperature + REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab + !$OMP THREADPRIVATE(tslab) + ! heat flux convergence due to Ekman + REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_ekman + !$OMP THREADPRIVATE(dt_ekman) + ! heat flux convergence due to horiz diffusion + REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_hdiff + !$OMP THREADPRIVATE(dt_hdiff) + ! heat flux convergence due to GM eddy advection + REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_gm + !$OMP THREADPRIVATE(dt_gm) + ! Heat Flux correction + REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: dt_qflux + !$OMP THREADPRIVATE(dt_qflux) + ! fraction of ocean covered by sea ice (sic / (oce+sic)) + REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic + !$OMP THREADPRIVATE(fsic) + ! temperature of the sea ice +!GG + !REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice_slab + !!$OMP THREADPRIVATE(tice_slab) + REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice_slab + !$OMP THREADPRIVATE(tice_slab) +!GG + ! sea ice thickness, in kg/m2 + REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice + !$OMP THREADPRIVATE(seaice) + ! net surface heat flux, weighted by open ocean fraction + ! slab_bils accumulated over cpl_pas timesteps + REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum + !$OMP THREADPRIVATE(bils_cum) + ! net heat flux into the ocean below the ice : conduction + solar radiation + REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg + !$OMP THREADPRIVATE(slab_bilg) + ! slab_bilg over cpl_pas timesteps + REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum + !$OMP THREADPRIVATE(bilg_cum) + ! wind stress saved over cpl_pas timesteps + REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: taux_cum + !$OMP THREADPRIVATE(taux_cum) + REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tauy_cum + !$OMP THREADPRIVATE(tauy_cum) + +!*********************************************************************************** +! Parameters (could be read in def file: move to slab_init) +!*********************************************************************************** +! snow and ice physical characteristics: + REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp + REAL, PARAMETER :: t_melt=273.15 ! melting ice temp + REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3 + REAL, PARAMETER :: ice_den=917. ! ice density + REAL, PARAMETER :: sea_den=1025. ! sea water density + REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice + REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow + REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice + REAL, PARAMETER :: sea_cap=3995. ! specific heat capacity, snow and ice + REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice + +! control of snow and ice cover & freeze / melt (heights converted to kg/m2) + REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm + REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean + REAL, PARAMETER :: ice_frac_min=0.001 + REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction + REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness + REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness + ! below ice_thin, priority is melt lateral / grow height + ! ice_thin is also height of new ice + REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness + ! above ice_thick, priority is melt height / grow lateral + REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice + REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height + +! albedo and radiation parameters + REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo + REAL, PARAMETER :: alb_sno_del=0.3 !max snow albedo = min + del + REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice + REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice + REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the + ! ice (no snow) + REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1) + +! horizontal transport + LOGICAL, PUBLIC, SAVE :: slab_hdiff + !$OMP THREADPRIVATE(slab_hdiff) + LOGICAL, PUBLIC, SAVE :: slab_gm + !$OMP THREADPRIVATE(slab_gm) + REAL, PRIVATE, SAVE :: coef_hdiff ! coefficient for horizontal diffusion + !$OMP THREADPRIVATE(coef_hdiff) + INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment + !$OMP THREADPRIVATE(slab_ekman) + !$OMP THREADPRIVATE(slab_cadj) + +!*********************************************************************************** + +CONTAINS +! +!*********************************************************************************** +! + SUBROUTINE ocean_slab_init(dtime, pctsrf_rst) + !, seaice_rst etc + + USE ioipsl_getin_p_mod, ONLY : getin_p + USE mod_phys_lmdz_transfert_para, ONLY : gather + USE slab_heat_transp_mod, ONLY : ini_slab_transp + + ! Input variables +!*********************************************************************************** + REAL, INTENT(IN) :: dtime +! Variables read from restart file + REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst + ! surface fractions from start file + +! Local variables +!************************************************************************************ + INTEGER :: error + REAL, DIMENSION(klon_glo) :: zmasq_glo + CHARACTER (len = 80) :: abort_message + CHARACTER (len = 20) :: modname = 'ocean_slab_intit' + +!*********************************************************************************** +! Define some parameters +!*********************************************************************************** +! Number of slab layers + nslay=2 + CALL getin_p('slab_layers',nslay) + print *,'number of slab layers : ',nslay +! Layer thickness + ALLOCATE(slabh(nslay), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation slabh' + CALL abort_physic(modname,abort_message,1) + ENDIF + slabh(1)=50. + CALL getin_p('slab_depth',slabh(1)) + IF (nslay.GT.1) THEN + slabh(2)=150. + END IF + +! cyang = 1/heat capacity of top layer (rho.c.H) + cyang=1/(slabh(1)*sea_den*sea_cap) + +! cpl_pas coupling period (update of tslab and ice fraction) +! pour un calcul a chaque pas de temps, cpl_pas=1 + cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour + CALL getin_p('cpl_pas',cpl_pas) + print *,'cpl_pas',cpl_pas + +! Horizontal diffusion + slab_hdiff=.FALSE. + CALL getin_p('slab_hdiff',slab_hdiff) + coef_hdiff=25000. + CALL getin_p('coef_hdiff',coef_hdiff) +! Ekman transport + slab_ekman=0 + CALL getin_p('slab_ekman',slab_ekman) +! GM eddy advection (2-layers only) + slab_gm=.FALSE. + CALL getin_p('slab_gm',slab_gm) + IF (slab_ekman.LT.2) THEN + slab_gm=.FALSE. + ENDIF +! Convective adjustment + IF (nslay.EQ.1) THEN + slab_cadj=0 + ELSE + slab_cadj=1 + END IF + CALL getin_p('slab_cadj',slab_cadj) + +!************************************************************************************ +! Allocate surface fraction read from restart file +!************************************************************************************ + ALLOCATE(fsic(klon), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation tmp_pctsrf_slab' + CALL abort_physic(modname,abort_message,1) + ENDIF + fsic(:)=0. + !zmasq = continent fraction + WHERE (1.-zmasq(:)>EPSFRA) + fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:)) + END WHERE + +!************************************************************************************ +! Allocate saved fields +!************************************************************************************ + ALLOCATE(tslab(klon,nslay), stat=error) + IF (error /= 0) CALL abort_physic & + (modname,'pb allocation tslab', 1) + + ALLOCATE(bils_cum(klon), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation slab_bils_cum' + CALL abort_physic(modname,abort_message,1) + ENDIF + bils_cum(:) = 0.0 + + IF (version_ocean=='sicINT') THEN ! interactive sea ice + ALLOCATE(slab_bilg(klon), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation slab_bilg' + CALL abort_physic(modname,abort_message,1) + ENDIF + slab_bilg(:) = 0.0 + ALLOCATE(bilg_cum(klon), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation slab_bilg_cum' + CALL abort_physic(modname,abort_message,1) + ENDIF + bilg_cum(:) = 0.0 +!GG ALLOCATE(tice(klon), stat = error) + ALLOCATE(tice_slab(klon), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation slab_tice' + CALL abort_physic(modname,abort_message,1) + ENDIF + ALLOCATE(seaice(klon), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation slab_seaice' + CALL abort_physic(modname,abort_message,1) + ENDIF + END IF + + IF (slab_hdiff) THEN !horizontal diffusion + ALLOCATE(dt_hdiff(klon,nslay), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation dt_hdiff' + CALL abort_physic(modname,abort_message,1) + ENDIF + dt_hdiff(:,:) = 0.0 + ENDIF + + ALLOCATE(dt_qflux(klon,nslay), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation dt_qflux' + CALL abort_physic(modname,abort_message,1) + ENDIF + dt_qflux(:,:) = 0.0 + + IF (slab_gm) THEN !GM advection + ALLOCATE(dt_gm(klon,nslay), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation dt_gm' + CALL abort_physic(modname,abort_message,1) + ENDIF + dt_gm(:,:) = 0.0 + ENDIF + + IF (slab_ekman.GT.0) THEN ! ekman transport + ALLOCATE(dt_ekman(klon,nslay), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation dt_ekman' + CALL abort_physic(modname,abort_message,1) + ENDIF + dt_ekman(:,:) = 0.0 + ALLOCATE(taux_cum(klon), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation taux_cum' + CALL abort_physic(modname,abort_message,1) + ENDIF + taux_cum(:) = 0.0 + ALLOCATE(tauy_cum(klon), stat = error) + IF (error /= 0) THEN + abort_message='Pb allocation tauy_cum' + CALL abort_physic(modname,abort_message,1) + ENDIF + tauy_cum(:) = 0.0 + ENDIF + +! Initialize transport + IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN + CALL gather(zmasq,zmasq_glo) +! Master thread/process only +!$OMP MASTER + IF (is_mpi_root) THEN + CALL ini_slab_transp(zmasq_glo) + END IF +!$OMP END MASTER + END IF + + END SUBROUTINE ocean_slab_init +! +!*********************************************************************************** +! + SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified) + +! this routine sends back the sea ice and ocean fraction to the main physics +! routine. Called only with interactive sea ice + +! Arguments +!************************************************************************************ + INTEGER, INTENT(IN) :: itime ! current timestep + INTEGER, INTENT(IN) :: jour ! day in year (not + REAL , INTENT(IN) :: dtime ! physics timestep (s) + REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg ! sub-surface fraction + LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is + ! modified at this time step + + pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:)) + pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:)) + is_modified=.TRUE. + + END SUBROUTINE ocean_slab_frac +! +!************************************************************************************ +! + SUBROUTINE ocean_slab_noice( & + itime, dtime, jour, knon, knindex, & + p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, & + AcoefH, AcoefQ, BcoefH, BcoefQ, & + AcoefU, AcoefV, BcoefU, BcoefV, & + ps, u1, v1, gustiness, tsurf_in, & + radsol, snow, & + qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & + tsurf_new, dflux_s, dflux_l, slab_bils) + + USE clesphys_mod_h + USE calcul_fluxs_mod + USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman1,slab_ekman2,slab_gmdiff + USE mod_phys_lmdz_para + + +! This routine +! (1) computes surface turbulent fluxes over points with some open ocean +! (2) reads additional Q-flux (everywhere) +! (3) computes horizontal transport (diffusion & Ekman) +! (4) updates slab temperature every cpl_pas ; creates new ice if needed. + +! Note : +! klon total number of points +! knon number of points with open ocean (varies with time) +! knindex gives position of the knon points within klon. +! In general, local saved variables have klon values +! variables exchanged with PBL module have knon. + +! Input arguments +!*********************************************************************************** + INTEGER, INTENT(IN) :: itime ! current timestep INTEGER, + INTEGER, INTENT(IN) :: jour ! day in year (for Q-Flux) + INTEGER, INTENT(IN) :: knon ! number of points + INTEGER, DIMENSION(klon), INTENT(IN) :: knindex + REAL, INTENT(IN) :: dtime ! timestep (s) + REAL, DIMENSION(klon), INTENT(IN) :: p1lay + REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm + ! drag coefficients + REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow + REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum ! near surface T, q + REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ + REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV + ! exchange coefficients for boundary layer scheme + REAL, DIMENSION(klon), INTENT(IN) :: ps ! surface pressure + REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! surface wind + REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! surface temperature + REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux + +! In/Output arguments +!************************************************************************************ + REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2 + +! Output arguments +!************************************************************************************ + REAL, DIMENSION(klon), INTENT(OUT) :: qsurf + REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat + REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 + REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! new surface tempearture + REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l + REAL, DIMENSION(klon), INTENT(OUT) :: slab_bils + +! Local variables +!************************************************************************************ + INTEGER :: i,ki,k + REAL :: t_cadj + ! for surface heat fluxes + REAL, DIMENSION(klon) :: cal, beta, dif_grnd + ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes + REAL, DIMENSION(klon) :: diff_sst, diff_siv + REAL, DIMENSION(klon,nslay) :: lmt_bils + ! for surface wind stress + REAL, DIMENSION(klon) :: u0, v0 + REAL, DIMENSION(klon) :: u1_lay, v1_lay + ! for new ice creation + REAL :: e_freeze, h_new, dfsic + ! horizontal diffusion and Ekman local vars + ! dimension = global domain (klon_glo) instead of // subdomains + REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo + ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop + REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp + REAL, DIMENSION(klon_glo,nslay) :: tslab_glo + REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo + +!**************************************************************************************** +! 1) Surface fluxes calculation +! +!**************************************************************************************** + !cal(:) = 0. ! infinite thermal inertia + !beta(:) = 1. ! wet surface + !dif_grnd(:) = 0. ! no diffusion into ground + ! EV: use calbeta + CALL calbeta(dtime, is_oce, knon, snow,qsurf, beta, cal, dif_grnd) + + + +! Suppose zero surface speed + u0(:)=0.0 + v0(:)=0.0 + u1_lay(:) = u1(:) - u0(:) + v1_lay(:) = v1(:) - v0(:) + +! Compute latent & sensible fluxes + CALL calcul_fluxs(knon, is_oce, dtime, & + tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, & + precip_rain, precip_snow, snow, qsurf, & + radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & + f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & + tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) + +! save total cumulated heat fluxes locally +! radiative + turbulent + melt of falling snow + slab_bils(:)=0. + DO i=1,knon + ki=knindex(i) + slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) & + -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki))) + bils_cum(ki)=bils_cum(ki)+slab_bils(ki) + END DO + +! Compute surface wind stress + CALL calcul_flux_wind(knon, dtime, & + u0, v0, u1, v1, gustiness, cdragm, & + AcoefU, AcoefV, BcoefU, BcoefV, & + p1lay, temp_air, & + flux_u1, flux_v1) + +! save cumulated wind stress + IF (slab_ekman.GT.0) THEN + DO i=1,knon + ki=knindex(i) + taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas + tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas + END DO + ENDIF + +!**************************************************************************************** +! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc +! +!**************************************************************************************** + CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) + ! lmt_bils and diff_sst,siv saved by limit_slab + ! qflux = total QFlux correction (in W/m2) + dt_qflux(:,1)=lmt_bils(:,1)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400. + IF (nslay.GT.1) THEN + dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay) + END IF + +!**************************************************************************************** +! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport) +! Bring to freezing temp and make sea ice if necessary +! +!***********************************************o***************************************** + tsurf_new=tsurf_in + IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction +! *********************************** +! Horizontal transport +! *********************************** + IF (slab_ekman.GT.0) THEN + ! copy wind stress to global var + CALL gather(taux_cum,taux_glo) + CALL gather(tauy_cum,tauy_glo) + END IF + + IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN + CALL gather(tslab,tslab_glo) + ! Compute horiz transport on one process only + IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus + IF (slab_hdiff) THEN + dt_hdiff_glo(:,:)=0. + END IF + IF (slab_ekman.GT.0) THEN + dt_ekman_glo(:,:)=0. + END IF + IF (slab_gm) THEN + dt_gm_glo(:,:)=0. + END IF + DO i=1,cpl_pas ! time splitting for numerical stability + IF (slab_ekman.GT.0) THEN + SELECT CASE (slab_ekman) + CASE (1) + CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp) + CASE (2) + CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm) + CASE DEFAULT + dt_ekman_tmp(:,:)=0. + END SELECT + dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:) + ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1 + DO k=1,nslay + dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den) + ENDDO + tslab_glo=tslab_glo+dt_ekman_tmp*dtime + IF (slab_gm) THEN ! Gent-McWilliams eddy advection + dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:) + ! convert dt from K.s-1.(kg.m-2) to K.s-1 + DO k=1,nslay + dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den) + END DO + tslab_glo=tslab_glo+dt_hdiff_tmp*dtime + END IF + ENDIF +! GM included in Ekman_2 +! IF (slab_gm) THEN ! Gent-McWilliams eddy advection +! CALL slab_gmdiff(tslab_glo,dt_hdiff_tmp) +! ! convert dt_gm from K.m.s-1 to K.s-1 +! DO k=1,nslay +! dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/slabh(k) +! END DO +! tslab_glo=tslab_glo+dt_hdiff_tmp*dtime +! dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:) +! END IF + IF (slab_hdiff) THEN ! horizontal diffusion + ! laplacian of slab T + CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp) + ! multiply by diff coef and normalize to 50m slab equivalent + dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh) + dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:) + tslab_glo=tslab_glo+dt_hdiff_tmp*dtime + END IF + END DO ! time splitting + IF (slab_hdiff) THEN + !dt_hdiff_glo saved in W/m2 + DO k=1,nslay + dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas + END DO + END IF + IF (slab_gm) THEN + !dt_hdiff_glo saved in W/m2 + dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas + END IF + IF (slab_ekman.GT.0) THEN + ! dt_ekman_glo saved in W/m2 + dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas + END IF + END IF ! master process +!$OMP BARRIER + ! Send new fields back to all processes + CALL Scatter(tslab_glo,tslab) + IF (slab_hdiff) THEN + CALL Scatter(dt_hdiff_glo,dt_hdiff) + END IF + IF (slab_gm) THEN + CALL Scatter(dt_gm_glo,dt_gm) + END IF + IF (slab_ekman.GT.0) THEN + CALL Scatter(dt_ekman_glo,dt_ekman) + ! clear wind stress + taux_cum(:)=0. + tauy_cum(:)=0. + END IF + ENDIF ! transport + +! *********************************** +! Other heat fluxes +! *********************************** + ! Add read QFlux + DO k=1,nslay + tslab(:,k)=tslab(:,k)+dt_qflux(:,k)*cyang*dtime*cpl_pas & + *slabh(1)/slabh(k) + END DO + ! Add cumulated surface fluxes + tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime + ! Convective adjustment if 2 layers + IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN + DO i=1,klon + IF (tslab(i,2).GT.tslab(i,1)) THEN + ! mean (mass-weighted) temperature + t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:)) + tslab(i,1)=t_cadj + tslab(i,2)=t_cadj + END IF + END DO + END IF +! *********************************** +! Update surface temperature and ice +! *********************************** + SELECT CASE(version_ocean) + CASE('sicNO') ! no sea ice even below freezing ! + DO i=1,knon + ki=knindex(i) + tsurf_new(i)=tslab(ki,1) + END DO + CASE('sicOBS') ! "realistic" case, for prescribed sea ice + ! tslab cannot be below freezing, or above it if there is sea ice + DO i=1,knon + ki=knindex(i) + IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN + tslab(ki,1)=t_freeze + END IF + tsurf_new(i)=tslab(ki,1) + END DO + CASE('sicINT') ! interactive sea ice + DO i=1,knon + ki=knindex(i) + IF (fsic(ki).LT.epsfra) THEN ! Free of ice + IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice + ! quantity of new ice formed + e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat + ! new ice + !GG tice(ki)=t_freeze + tice_slab(ki)=t_freeze + fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin) + IF (fsic(ki).GT.ice_frac_min) THEN + seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max) + tslab(ki,1)=t_freeze + ELSE + fsic(ki)=0. + END IF + tsurf_new(i)=t_freeze + ELSE + tsurf_new(i)=tslab(ki,1) + END IF + ELSE ! ice present + tsurf_new(i)=t_freeze + IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice + ! quantity of new ice formed over open ocean + e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) & + /(ice_lat+ice_cap/2.*(t_freeze-tice_slab(ki))) + !GG /(ice_lat+ice_cap/2.*(t_freeze-tice(ki))) + ! new ice height and fraction + h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new + dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new) + h_new=MIN(e_freeze/dfsic,h_ice_max) + ! update tslab to freezing over open ocean only + tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki)) + ! update sea ice + seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) & + /(dfsic+fsic(ki)) + fsic(ki)=fsic(ki)+dfsic + ! update snow? + END IF ! tslab below freezing + END IF ! sea ice present + END DO + END SELECT + bils_cum(:)=0.0! clear cumulated fluxes + END IF ! coupling time + END SUBROUTINE ocean_slab_noice +! +!**************************************************************************************** + + SUBROUTINE ocean_slab_ice( & + itime, dtime, jour, knon, knindex, & + tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, & + AcoefH, AcoefQ, BcoefH, BcoefQ, & + AcoefU, AcoefV, BcoefU, BcoefV, & + ps, u1, v1, gustiness, & + radsol, snow, qsurf, qsol, agesno, & + alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & + tsurf_new, dflux_s, dflux_l, swnet) + + USE clesphys_mod_h + USE yomcst_mod_h +USE calcul_fluxs_mod + + + +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: itime, jour, knon + INTEGER, DIMENSION(klon), INTENT(IN) :: knindex + REAL, INTENT(IN) :: dtime + REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in + REAL, DIMENSION(klon), INTENT(IN) :: p1lay + REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm + REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow + REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum + REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ + REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV + REAL, DIMENSION(klon), INTENT(IN) :: ps + REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness + REAL, DIMENSION(klon), INTENT(IN) :: swnet + +! In/Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol + REAL, DIMENSION(klon), INTENT(INOUT) :: agesno + REAL, DIMENSION(klon), INTENT(INOUT) :: radsol + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(OUT) :: qsurf + REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval + REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval + REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat + REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 + REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new + REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l + +! Local variables +!**************************************************************************************** + INTEGER :: i,ki + REAL, DIMENSION(klon) :: cal, beta, dif_grnd + REAL, DIMENSION(klon) :: u0, v0 + REAL, DIMENSION(klon) :: u1_lay, v1_lay + ! intermediate heat fluxes: + REAL :: f_cond, f_swpen + ! for snow/ice albedo: + REAL :: alb_snow, alb_ice, alb_pond + REAL :: frac_snow, frac_ice, frac_pond + ! for ice melt / freeze + REAL :: e_melt, snow_evap, h_test + ! dhsic, dfsic change in ice mass, fraction. + REAL :: dhsic, dfsic, frac_mf + +!**************************************************************************************** +! 1) Flux calculation +!**************************************************************************************** +! Suppose zero surface speed + u0(:)=0.0 + v0(:)=0.0 + u1_lay(:) = u1(:) - u0(:) + v1_lay(:) = v1(:) - v0(:) + +! set beta, cal, compute conduction fluxes inside ice/snow + slab_bilg(:)=0. + !dif_grnd(:)=0. + !beta(:) = 1. + ! EV: use calbeta to calculate beta and then recalculate properly cal + CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, cal, dif_grnd) + + + DO i=1,knon + ki=knindex(i) + IF (snow(i).GT.snow_min) THEN + ! snow-layer heat capacity + cal(i)=2.*RCPD/(snow(i)*ice_cap) + ! snow conductive flux + f_cond=sno_cond*(tice_slab(ki)-tsurf_in(i))/snow(i) + !GG f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i) + ! all shortwave flux absorbed + f_swpen=0. + ! bottom flux (ice conduction) + slab_bilg(ki)=ice_cond*(tice_slab(ki)-t_freeze)/seaice(ki) + !GG slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki) + ! update ice temperature + !GG tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) & + tice_slab(ki)=tice_slab(ki)-2./ice_cap/(snow(i)+seaice(ki)) & + *(slab_bilg(ki)+f_cond)*dtime + ELSE ! bare ice + ! ice-layer heat capacity + cal(i)=2.*RCPD/(seaice(ki)*ice_cap) + ! conductive flux + !GG f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki) + f_cond=ice_cond*(t_freeze-tice_slab(ki))/seaice(ki) + ! penetrative shortwave flux... + f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den) + slab_bilg(ki)=f_swpen-f_cond + END IF + radsol(i)=radsol(i)+f_cond-f_swpen + END DO + ! weight fluxes to ocean by sea ice fraction + slab_bilg(:)=slab_bilg(:)*fsic(:) + +! calcul_fluxs (sens, lat etc) + CALL calcul_fluxs(knon, is_sic, dtime, & + tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, & + precip_rain, precip_snow, snow, qsurf, & + radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & + f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & + tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) + DO i=1,knon + !GG IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i) + IF (snow(i).LT.snow_min) tice_slab(knindex(i))=tsurf_new(i) + END DO + +! calcul_flux_wind + CALL calcul_flux_wind(knon, dtime, & + u0, v0, u1, v1, gustiness, cdragm, & + AcoefU, AcoefV, BcoefU, BcoefV, & + p1lay, temp_air, & + flux_u1, flux_v1) + +!**************************************************************************************** +! 2) Update snow and ice surface +!**************************************************************************************** +! snow precip + DO i=1,knon + ki=knindex(i) + IF (precip_snow(i) > 0.) THEN + snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki))) + END IF +! snow and ice sublimation + IF (evap(i) > 0.) THEN + snow_evap = MIN (snow(i) / dtime, evap(i)) + snow(i) = snow(i) - snow_evap * dtime + snow(i) = MAX(0.0, snow(i)) + seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime) + ENDIF +! Melt / Freeze snow from above if Tsurf>0 + IF (tsurf_new(i).GT.t_melt) THEN + ! energy available for melting snow (in kg of melted snow /m2) + e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. & + /(ice_lat+ice_cap/2.*(t_melt-tice_slab(ki))),0.0),snow(i)) + !GG /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i)) + ! remove snow + IF (snow(i).GT.e_melt) THEN + snow(i)=snow(i)-e_melt + tsurf_new(i)=t_melt + ELSE ! all snow is melted + ! add remaining heat flux to ice + e_melt=e_melt-snow(i) + tice_slab(ki)=tice_slab(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki)) + !GG tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki)) + tsurf_new(i)=tice_slab(ki) + !GG tsurf_new(i)=tice(ki) + END IF + END IF +! melt ice from above if Tice>0 + !GG IF (tice(ki).GT.t_melt) THEN + IF (tice_slab(ki).GT.t_melt) THEN + ! quantity of ice melted (kg/m2) + !GG e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. & + e_melt=MAX(seaice(ki)*(tice_slab(ki)-t_melt)*ice_cap/2. & + /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0) + ! melt from above, height only + dhsic=MIN(seaice(ki)-h_ice_min,e_melt) + e_melt=e_melt-dhsic + IF (e_melt.GT.0) THEN + ! lateral melt if ice too thin + dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki)) + ! if all melted add remaining heat to ocean + e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min) + slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime + ! update height and fraction + fsic(ki)=fsic(ki)-dfsic + END IF + seaice(ki)=seaice(ki)-dhsic + ! surface temperature at melting point + !GG tice(ki)=t_melt + tice_slab(ki)=t_melt + tsurf_new(i)=t_melt + END IF + ! convert snow to ice if below floating line + h_test=(seaice(ki)+snow(i))*ice_den-seaice(ki)*sea_den + IF (h_test.GT.0.) THEN !snow under water + ! extra snow converted to ice (with added frozen sea water) + dhsic=h_test/(sea_den-ice_den+sno_den) + seaice(ki)=seaice(ki)+dhsic + snow(i)=snow(i)-dhsic*sno_den/ice_den + ! available energy (freeze sea water + bring to tice_slab) + e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ & + ice_cap/2.*(t_freeze-tice_slab(ki))) + !GG ice_cap/2.*(t_freeze-tice(ki))) + ! update ice temperature + !GG tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki)) + tice_slab(ki)=tice_slab(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki)) + END IF + END DO + +! New albedo + DO i=1,knon + ki=knindex(i) + ! snow albedo: update snow age + IF (snow(i).GT.0.0001) THEN + agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)& + * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.) + ELSE + agesno(i)=0.0 + END IF + ! snow albedo + alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.) + ! ice albedo (varies with ice tkickness and temp) + alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1) + !GG IF (tice(ki).GT.t_freeze-0.01) THEN + IF (tice_slab(ki).GT.t_freeze-0.01) THEN + alb_ice=MIN(alb_ice,alb_ice_wet) + ELSE + alb_ice=MIN(alb_ice,alb_ice_dry) + END IF + ! pond albedo + alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice_slab(ki)-t_melt,-2.0))) + !GG alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) + ! pond fraction + frac_pond=0.2*(2.0+MIN(0.0,MAX(tice_slab(ki)-t_melt,-2.0))) + !GG frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) + ! snow fraction + frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min)) + ! ice fraction + frac_ice=MAX(0.0,1.-frac_pond-frac_snow) + ! total albedo + alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond + END DO + alb2_new(:) = alb1_new(:) + +!**************************************************************************************** +! 3) Recalculate new ocean temperature (add fluxes below ice) +! Melt / freeze from below +!***********************************************o***************************************** + !cumul fluxes + bilg_cum(:)=bilg_cum(:)+slab_bilg(:) + IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction + ! Add cumulated surface fluxes + tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime + DO i=1,knon + ki=knindex(i) + ! split lateral/top melt-freeze + frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin))) + IF (tslab(ki,1).LE.t_freeze) THEN + ! ****** Form new ice from below ******* + ! quantity of new ice + e_melt=(t_freeze-tslab(ki,1))/cyang & + /(ice_lat+ice_cap/2.*(t_freeze-tice_slab(ki))) + !GG /(ice_lat+ice_cap/2.*(t_freeze-tice(ki))) + ! first increase height to h_thin + dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki))) + seaice(ki)=dhsic+seaice(ki) + e_melt=e_melt-fsic(ki)*dhsic + IF (e_melt.GT.0.) THEN + ! frac_mf fraction used for lateral increase + dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki)) + fsic(ki)=fsic(ki)+dfsic + e_melt=e_melt-dfsic*seaice(ki) + ! rest used to increase height + seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki)) + END IF + tslab(ki,1)=t_freeze + ELSE ! slab temperature above freezing + ! ****** melt ice from below ******* + ! quantity of melted ice + e_melt=(tslab(ki,1)-t_freeze)/cyang & + /(ice_lat+ice_cap/2.*(tice_slab(ki)-t_freeze)) + !GG /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze)) + ! first decrease height to h_thick + dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki))) + seaice(ki)=seaice(ki)-dhsic + e_melt=e_melt-fsic(ki)*dhsic + IF (e_melt.GT.0) THEN + ! frac_mf fraction used for height decrease + dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki))) + seaice(ki)=seaice(ki)-dhsic + e_melt=e_melt-fsic(ki)*dhsic + ! rest used to decrease fraction (up to 0!) + dfsic=MIN(fsic(ki),e_melt/seaice(ki)) + ! keep remaining in ocean + e_melt=e_melt-dfsic*seaice(ki) + END IF + tslab(ki,1)=t_freeze+e_melt*ice_lat*cyang + fsic(ki)=fsic(ki)-dfsic + END IF + END DO + bilg_cum(:)=0. + END IF ! coupling time + + !tests ice fraction + WHERE (fsic.LT.ice_frac_min) + tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang + tice_slab=t_melt + !GG tice=t_melt + fsic=0. + seaice=0. + END WHERE + + END SUBROUTINE ocean_slab_ice +! +!**************************************************************************************** +! + SUBROUTINE ocean_slab_final + +!**************************************************************************************** +! Deallocate module variables +!**************************************************************************************** + IF (ALLOCATED(tslab)) DEALLOCATE(tslab) + IF (ALLOCATED(fsic)) DEALLOCATE(fsic) + IF (ALLOCATED(tice_slab)) DEALLOCATE(tice_slab) + IF (ALLOCATED(seaice)) DEALLOCATE(seaice) + IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg) + IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum) + IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum) + IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum) + IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum) + IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman) + IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff) + IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm) + IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux) + + END SUBROUTINE ocean_slab_final + +END MODULE ocean_slab_mod Index: LMDZ6/trunk/libf/phylmdiso/soil.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/soil.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/soil.f90 (revision 5986) @@ -0,0 +1,352 @@ +! +! $Header$ +! +SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, & + lon, lat, ptsoil, pcapcal, pfluxgrd) + + USE yomcst_mod_h + USE dimphy + USE mod_phys_lmdz_para + USE indice_sol_mod + USE print_control_mod, ONLY: lunout + USE dimsoil_mod_h, ONLY: nsoilmx + USE comsoil_mod_h + IMPLICIT NONE + +!======================================================================= +! +! Auteur: Frederic Hourdin 30/01/92 +! ------- +! +! Object: Computation of : the soil temperature evolution +! ------- the surfacic heat capacity "Capcal" +! the surface conduction flux pcapcal +! +! Update: 2021/07 : soil thermal inertia, formerly a constant value, +! ------ can also be now a function of soil moisture (F Cheruy's idea) +! depending on iflag_inertie, read from physiq.def via conf_phys_m.F90 +! ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy) +! +! Method: Implicit time integration +! ------- +! Consecutive ground temperatures are related by: +! T(k+1) = C(k) + D(k)*T(k) (*) +! The coefficients C and D are computed at the t-dt time-step. +! Routine structure: +! 1) C and D coefficients are computed from the old temperature +! 2) new temperatures are computed using (*) +! 3) C and D coefficients are computed from the new temperature +! profile for the t+dt time-step +! 4) the coefficients A and B are computed where the diffusive +! fluxes at the t+dt time-step is given by +! Fdiff = A + B Ts(t+dt) +! or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt +! with F0 = A + B (Ts(t)) +! Capcal = B*dt +! +! Interface: +! ---------- +! +! Arguments: +! ---------- +! ptimestep physical timestep (s) +! indice sub-surface index +! snow(klon) snow +! ptsrf(klon) surface temperature at time-step t (K) +! qsol(klon) soil moisture (kg/m2 or mm) +! lon(klon) longitude in radian +! lat(klon) latitude in radian +! ptsoil(klon,nsoilmx) temperature inside the ground (K) +! pcapcal(klon) surfacic specific heat (W*m-2*s*K-1) +! pfluxgrd(klon) surface diffusive flux from ground (Wm-2) +! +!======================================================================= +! Arguments +! --------- + REAL, INTENT(IN) :: ptimestep + INTEGER, INTENT(IN) :: indice, knon !, knindex + REAL, DIMENSION(klon), INTENT(IN) :: snow + REAL, DIMENSION(klon), INTENT(IN) :: ptsrf + REAL, DIMENSION(klon), INTENT(IN) :: qsol + REAL, DIMENSION(klon), INTENT(IN) :: lon + REAL, DIMENSION(klon), INTENT(IN) :: lat + + REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil + REAL, DIMENSION(klon), INTENT(OUT) :: pcapcal + REAL, DIMENSION(klon), INTENT(OUT) :: pfluxgrd + +!----------------------------------------------------------------------- +! Local variables +! --------------- + INTEGER :: ig, jk, ierr + REAL :: min_period,dalph_soil + REAL, DIMENSION(nsoilmx) :: zdz2 + REAL :: z1s + REAL, DIMENSION(klon) :: ztherm_i + REAL, DIMENSION(klon,nsoilmx,nbsrf) :: C_coef, D_coef + +! Local saved variables +! --------------------- + REAL, SAVE :: lambda +!$OMP THREADPRIVATE(lambda) + REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2 +!$OMP THREADPRIVATE(dz1,dz2) + LOGICAL, SAVE :: firstcall=.TRUE. +!$OMP THREADPRIVATE(firstcall) + +!----------------------------------------------------------------------- +! Depthts: +! -------- + REAL fz,rk,fz1,rk1,rk2 + fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) + + +!----------------------------------------------------------------------- +! Calculation of some constants +! NB! These constants do not depend on the sub-surfaces +!----------------------------------------------------------------------- + + IF (firstcall) THEN +!----------------------------------------------------------------------- +! ground levels +! grnd=z/l where l is the skin depth of the diurnal cycle: +!----------------------------------------------------------------------- + + min_period=1800. ! en secondes + dalph_soil=2. ! rapport entre les epaisseurs de 2 couches succ. +!$OMP MASTER + IF (is_mpi_root) THEN + OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr) + IF (ierr == 0) THEN ! Read file only if it exists + READ(99,*) min_period + READ(99,*) dalph_soil + WRITE(lunout,*)'Discretization for the soil model' + WRITE(lunout,*)'First level e-folding depth',min_period, & + ' dalph',dalph_soil + CLOSE(99) + END IF + ENDIF +!$OMP END MASTER + CALL bcast(min_period) + CALL bcast(dalph_soil) + +! la premiere couche represente un dixieme de cycle diurne + fz1=SQRT(min_period/3.14) + + DO jk=1,nsoilmx + rk1=jk + rk2=jk-1 + dz2(jk)=fz(rk1)-fz(rk2) + ENDDO + DO jk=1,nsoilmx-1 + rk1=jk+.5 + rk2=jk-.5 + dz1(jk)=1./(fz(rk1)-fz(rk2)) + ENDDO + lambda=fz(.5)*dz1(1) + WRITE(lunout,*)'full layers, intermediate layers (seconds)' + DO jk=1,nsoilmx + rk=jk + rk1=jk+.5 + rk2=jk-.5 + WRITE(lunout,*)'fz=', & + fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14 + ENDDO + + firstcall =.FALSE. + END IF + + +!----------------------------------------------------------------------- +! Calcul de l'inertie thermique a partir de la variable rnat. +! on initialise a inertie_sic meme au-dessus d'un point de mer au cas +! ou le point de mer devienne point de glace au pas suivant +! on corrige si on a un point de terre avec ou sans glace +! +! iophys can be used to write the ztherm_i variable in a phys.nc file +! and check the results; to do so, add "CALL iophys_ini" in physiq_mod +! and add knindex to the list of inputs in all the calls to soil.F90 +! (and to soil.F90 itself !) +!----------------------------------------------------------------------- + + IF (indice == is_sic) THEN + DO ig = 1, knon + ztherm_i(ig) = inertie_sic + ENDDO + IF (iflag_sic == 0) THEN + DO ig = 1, knon + IF (snow(ig) > 0.0) ztherm_i(ig) = inertie_sno + ENDDO +! Otherwise sea-ice keeps the same inertia, even when covered by snow + ENDIF +! CALL iophys_ecrit_index('ztherm_sic', 1, 'ztherm_sic', 'USI', & +! knon, knindex, ztherm_i) + ELSE IF (indice == is_lic) THEN + DO ig = 1, knon + ztherm_i(ig) = inertie_lic + IF (snow(ig) > 0.0) ztherm_i(ig) = inertie_sno + ENDDO +! CALL iophys_ecrit_index('ztherm_lic', 1, 'ztherm_lic', 'USI', & +! knon, knindex, ztherm_i) + ELSE IF (indice == is_ter) THEN + ! + ! La relation entre l'inertie thermique du sol et qsol change d'apres + ! iflag_inertie, defini dans physiq.def, et appele via comsoil.h + ! + DO ig = 1, knon + ! iflag_inertie=0 correspond au cas inertie=constant, comme avant + IF (iflag_inertie==0) THEN + ztherm_i(ig) = inertie_sol + ELSE IF (iflag_inertie == 1) THEN + ! I = a_qsol * qsol + b modele lineaire deduit d'une + ! regression lineaire I = a_mrsos * mrsos + b obtenue sur + ! sorties MO d'une simulation LMDZOR(CMIP6) sur l'annee 2000 + ! sur tous les points avec frac_snow=0 + ! Difference entre qsol et mrsos prise en compte par un + ! facteur d'echelle sur le coefficient directeur de regression: + ! fact = 35./150. = mrsos_max/qsol_max + ! et a_qsol = a_mrsos * fact (car a = dI/dHumidite) + ztherm_i(ig) = 30.0 *35.0/150.0 *qsol(ig) +770.0 + ! AS : pour qsol entre 0 - 150, on a I entre 770 - 1820 + ELSE IF (iflag_inertie == 2) THEN + ! deux regressions lineaires, sur les memes sorties, + ! distinguant le type de sol : sable ou autre (limons/argile) + ! Implementation simple : regression type "sable" seulement pour + ! Sahara, defini par une "boite" lat/lon (NB : en radians !! ) + IF (lon(ig)>-0.35 .AND. lon(ig)<0.70 .AND. lat(ig)>0.17 .AND. lat(ig)<0.52) THEN + ! Valeurs theoriquement entre 728 et 2373 ; qsol valeurs basses + ztherm_i(ig) = 47. *35.0/150.0 *qsol(ig) +728. ! boite type "sable" pour Sahara + ELSE + ! Valeurs theoriquement entre 550 et 1940 ; qsol valeurs moyennes et hautes + ztherm_i(ig) = 41. *35.0/150.0 *qsol(ig) +505. + ENDIF + ELSE IF (iflag_inertie == 3) THEN + ! AS : idee a tester : + ! si la relation doit etre une droite, + ! definissons-la en fonction des valeurs min et max de qsol (0:150), + ! et de l'inertie (900 : 2000 ou 2400 ; choix ici: 2000) + ! I = I_min + qsol * (I_max - I_min)/(qsol_max - qsol_min) + ztherm_i(ig) = 900. + qsol(ig) * (2000. - 900.)/150. + ELSE + WRITE (lunout,*) "Le choix iflag_inertie = ",iflag_inertie," n'est pas defini. Veuillez choisir un entier entre 0 et 3" + ENDIF + ! + ! Fin de l'introduction de la relation entre l'inertie thermique du sol et qsol + !------------------------------------------- + !AS : donc le moindre flocon de neige sur un point de grid + ! fait que l'inertie du point passe a la valeur pour neige ! + IF (snow(ig) > 0.0) ztherm_i(ig) = inertie_sno + + ENDDO +! CALL iophys_ecrit_index('ztherm_ter', 1, 'ztherm_ter', 'USI', & +! knon, knindex, ztherm_i) + ELSE IF (indice == is_oce) THEN + DO ig = 1, knon +! This is just in case, but SST should be used by the model anyway + ztherm_i(ig) = inertie_sic + ENDDO +! CALL iophys_ecrit_index('ztherm_oce', 1, 'ztherm_oce', 'USI', & +! knon, knindex, ztherm_i) + ELSE + WRITE(lunout,*) "valeur d indice non prevue", indice + call abort_physic("soil", "", 1) + ENDIF + + +!----------------------------------------------------------------------- +! 1) +! Calculation of Cgrf and Dgrd coefficients using soil temperature from +! previous time step. +! +! These variables are recalculated on the local compressed grid instead +! of saved in restart file. +!----------------------------------------------------------------------- + DO jk=1,nsoilmx + zdz2(jk)=dz2(jk)/ptimestep + ENDDO + + DO ig=1,knon + z1s = zdz2(nsoilmx)+dz1(nsoilmx-1) + C_coef(ig,nsoilmx-1,indice)= & + zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s + D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s + ENDDO + + DO jk=nsoilmx-1,2,-1 + DO ig=1,knon + z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) & + *(1.-D_coef(ig,jk,indice))) + C_coef(ig,jk-1,indice)= & + (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s + D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s + ENDDO + ENDDO + +!----------------------------------------------------------------------- +! 2) +! Computation of the soil temperatures using the Cgrd and Dgrd +! coefficient computed above +! +!----------------------------------------------------------------------- + +! Surface temperature + DO ig=1,knon + ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/ & + (lambda*(1.-D_coef(ig,1,indice))+1.) + ENDDO + +! Other temperatures + DO jk=1,nsoilmx-1 + DO ig=1,knon + ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) & + *ptsoil(ig,jk) + ENDDO + ENDDO + + IF (indice == is_sic) THEN + DO ig = 1 , knon + ptsoil(ig,nsoilmx) = RTT - 1.8 + END DO + ENDIF + +!----------------------------------------------------------------------- +! 3) +! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil +! temperature +!----------------------------------------------------------------------- + DO ig=1,knon + z1s = zdz2(nsoilmx)+dz1(nsoilmx-1) + C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s + D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s + ENDDO + + DO jk=nsoilmx-1,2,-1 + DO ig=1,knon + z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) & + *(1.-D_coef(ig,jk,indice))) + C_coef(ig,jk-1,indice) = & + (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s + D_coef(ig,jk-1,indice) = dz1(jk-1)*z1s + ENDDO + ENDDO + +!----------------------------------------------------------------------- +! 4) +! Computation of the surface diffusive flux from ground and +! calorific capacity of the ground +!----------------------------------------------------------------------- + DO ig=1,knon + pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* & + (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1)) + pcapcal(ig) = ztherm_i(ig)* & + (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1)) + z1s = lambda*(1.-D_coef(ig,1,indice))+1. + pcapcal(ig) = pcapcal(ig)/z1s + pfluxgrd(ig) = pfluxgrd(ig) & + + pcapcal(ig) * (ptsoil(ig,1) * z1s & + - lambda * C_coef(ig,1,indice) & + - ptsrf(ig)) & + /ptimestep + ENDDO + +END SUBROUTINE soil Index: LMDZ6/trunk/libf/phylmdiso/surf_seaice_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/surf_seaice_mod.F90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/surf_seaice_mod.F90 (revision 5986) @@ -0,0 +1,235 @@ +! +! $Id: surf_seaice_mod.F90 5662 2025-05-20 14:24:41Z fairhead $ +! +MODULE surf_seaice_mod + + IMPLICIT NONE + +CONTAINS +! +!**************************************************************************************** +! + SUBROUTINE surf_seaice( & + rlon, rlat, swnet, lwnet, alb1, fder, & + itime, dtime, jour, knon, knindex, & + lafin, & + tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, & + AcoefH, AcoefQ, BcoefH, BcoefQ, & + AcoefU, AcoefV, BcoefU, BcoefV, & + ps, u1, v1, gustiness, pctsrf, & + snow, qsurf, qsol, agesno, tsoil, & + z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, & + tsurf_new, dflux_s, dflux_l, & +!GG flux_u1, flux_v1) + flux_u1, flux_v1, hice, tice,bilg_cumul, & + fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, & + dtice_melt, dtice_snow2sic & +!GG +#ifdef ISO + & ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, & + & xtsnow,xtsol,xtevap,Rland_ice & +#endif + & ) + + USE dimphy + USE surface_data + USE ocean_forced_mod, ONLY : ocean_forced_ice + USE ocean_cpl_mod, ONLY : ocean_cpl_ice + USE ocean_slab_mod, ONLY : ocean_slab_ice + USE indice_sol_mod +#ifdef ISO + USE infotrac_phy, ONLY : ntiso,niso +#endif + USE clesphys_mod_h + USE yomcst_mod_h +USE dimsoil_mod_h, ONLY: nsoilmx + +! +! This subroutine will make a call to ocean_XXX_ice according to the ocean mode (force, +! slab or couple). The calculation of rugosity for the sea-ice surface is also done +! in here because it is the same calculation for the different modes of ocean. +! + + + ! for rd and retv + +! Input arguments +!**************************************************************************************** + INTEGER, INTENT(IN) :: itime, jour, knon + INTEGER, DIMENSION(klon), INTENT(IN) :: knindex + LOGICAL, INTENT(IN) :: lafin + REAL, INTENT(IN) :: dtime + REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat + REAL, DIMENSION(klon), INTENT(IN) :: swnet ! net shortwave radiation at surface + REAL, DIMENSION(klon), INTENT(IN) :: lwnet ! net longwave radiation at surface + REAL, DIMENSION(klon), INTENT(IN) :: alb1 ! albedo in visible SW interval + REAL, DIMENSION(klon), INTENT(IN) :: fder + REAL, DIMENSION(klon), INTENT(IN) :: tsurf + REAL, DIMENSION(klon), INTENT(IN) :: p1lay + REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm + REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow + REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum + REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ + REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV + REAL, DIMENSION(klon), INTENT(IN) :: ps + REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness + REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf +#ifdef ISO + REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow + REAL, DIMENSION(klon), INTENT(IN) :: xtspechum + REAL, DIMENSION(niso,klon), INTENT(IN) :: Roce + REAL, DIMENSION(niso,klon), INTENT(IN) :: Rland_ice +#endif + +! In/Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsurf, qsol + REAL, DIMENSION(klon), INTENT(INOUT) :: agesno + REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil +#ifdef ISO + REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsnow + REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsol +#endif + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon), INTENT(OUT) :: z0m, z0h +!albedo SB >>> +! REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval +! REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval + REAL, DIMENSION(6), INTENT(IN) :: SFRWL + REAL, DIMENSION(klon,nsw), INTENT(OUT) :: alb_dir_new,alb_dif_new +!albedo SB <<< + REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat + REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new + REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l + REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 +!GG + REAL, DIMENSION(klon), INTENT(INOUT) :: hice, tice, bilg_cumul + REAL, DIMENSION(klon), INTENT(INOUT) :: fcds,fcdi, dh_basal_growth,dh_basal_melt + REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt, dh_snow2sic, dtice_melt, dtice_snow2sic +!GG +#ifdef ISO + REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap +#endif + +! Local arguments +!**************************************************************************************** + REAL, DIMENSION(klon) :: radsol +#ifdef ISO +#ifdef ISOVERIF + INTEGER :: j +#endif +#endif + +!albedo SB >>> + REAL, DIMENSION(klon) :: alb1_new,alb2_new +!albedo SB <<< + + real rhoa(knon) ! density of moist air (kg / m3) + +! End definitions +!**************************************************************************************** + + +!**************************************************************************************** +! Calculate total net radiance at surface +! +!**************************************************************************************** + radsol(:) = 0.0 + radsol(1:knon) = swnet(1:knon) + lwnet(1:knon) + + rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon))) + +!**************************************************************************************** +! Switch according to type of ocean (couple, slab or forced) +! +!**************************************************************************************** + IF (type_ocean == 'couple') THEN + + CALL ocean_cpl_ice( & + rlon, rlat, swnet, lwnet, alb1, & + fder, & + itime, dtime, knon, knindex, & + lafin,& + p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum,& + AcoefH, AcoefQ, BcoefH, BcoefQ, & + AcoefU, AcoefV, BcoefU, BcoefV, & + ps, u1, v1, gustiness, pctsrf, & + radsol, snow, qsurf, & + alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & + tsurf_new, dflux_s, dflux_l, rhoa) + + ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN + CALL ocean_slab_ice( & + itime, dtime, jour, knon, knindex, & + tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum,& + AcoefH, AcoefQ, BcoefH, BcoefQ, & + AcoefU, AcoefV, BcoefU, BcoefV, & + ps, u1, v1, gustiness, & + radsol, snow, qsurf, qsol, agesno, & + alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & + tsurf_new, dflux_s, dflux_l, swnet) + + ELSE ! type_ocean=force or slab +sicOBS or sicNO + CALL ocean_forced_ice( & + itime, dtime, jour, knon, knindex, & + tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum,& + AcoefH, AcoefQ, BcoefH, BcoefQ, & + AcoefU, AcoefV, BcoefU, BcoefV, & +!GG ps, u1, v1, gustiness, & + ps, u1, v1, gustiness,pctsrf, & +!GG + radsol, snow, qsol, agesno, tsoil, & + qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & +!GG tsurf_new, dflux_s, dflux_l, rhoa) + tsurf_new, dflux_s, dflux_l,rhoa,swnet,hice, tice, bilg_cumul, & + fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, & + dtice_melt, dtice_snow2sic & +!GG +#ifdef ISO + ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, & + xtsnow, xtsol,xtevap,Rland_ice & +#endif + ) + + END IF + +!**************************************************************************************** +! Calculate rugosity +! +!**************************************************************************************** + + z0m=z0m_seaice + z0h = z0h_seaice + +!albedo SB >>> + select case(NSW) + case(2) + alb_dir_new(1:knon,1)=alb1_new(1:knon) + alb_dir_new(1:knon,2)=alb2_new(1:knon) + case(4) + alb_dir_new(1:knon,1)=alb1_new(1:knon) + alb_dir_new(1:knon,2)=alb2_new(1:knon) + alb_dir_new(1:knon,3)=alb2_new(1:knon) + alb_dir_new(1:knon,4)=alb2_new(1:knon) + case(6) + alb_dir_new(1:knon,1)=alb1_new(1:knon) + alb_dir_new(1:knon,2)=alb1_new(1:knon) + alb_dir_new(1:knon,3)=alb1_new(1:knon) + alb_dir_new(1:knon,4)=alb2_new(1:knon) + alb_dir_new(1:knon,5)=alb2_new(1:knon) + alb_dir_new(1:knon,6)=alb2_new(1:knon) + end select +alb_dif_new=alb_dir_new +!albedo SB <<< + + + + + END SUBROUTINE surf_seaice +! +!**************************************************************************************** +! +END MODULE surf_seaice_mod + Index: LMDZ6/trunk/libf/phylmdiso/ustarhb.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/ustarhb.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/ustarhb.f90 (revision 5986) @@ -0,0 +1,43 @@ + +! $Header$ + +SUBROUTINE ustarhb(klon, klev, knon, u, v, cd_m, ustar) + IMPLICIT NONE + ! ====================================================================== + ! Laurent Li (LMD/CNRS), le 30 septembre 1998 + ! Couche limite non-locale. Adaptation du code du CCM3. + ! Code non teste, donc a ne pas utiliser. + ! ====================================================================== + ! Nonlocal scheme that determines eddy diffusivities based on a + ! diagnosed boundary layer height and a turbulent velocity scale. + ! Also countergradient effects for heat and moisture are included. + + ! For more information, see Holtslag, A.A.M., and B.A. Boville, 1993: + ! Local versus nonlocal boundary-layer diffusion in a global climate + ! model. J. of Climate, vol. 6, 1825-1842. + ! ====================================================================== + + ! Arguments: + + INTEGER, INTENT(IN) :: klon, klev, knon ! nombre de points a calculer + REAL, DIMENSION(klon, klev), INTENT(IN) :: u,v ! vent horizontal (m/s) + REAL, DIMENSION(klon), INTENT(IN) :: cd_m ! coefficient de friction au sol pour vitesse + REAL, DIMENSION(klon), INTENT(OUT) :: ustar + + INTEGER :: i, k + REAL :: zxt, zxq, zxu, zxv, zxmod, taux, tauy + REAL :: zx_alf1, zx_alf2 ! parametres pour extrapolation + + DO i = 1, knon + zx_alf1 = 1.0 + zx_alf2 = 1.0 - zx_alf1 + zxu = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2 + zxv = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2 + zxmod = 1.0 + sqrt(zxu**2+zxv**2) + taux = zxu*zxmod*cd_m(i) + tauy = zxv*zxmod*cd_m(i) + ustar(i) = sqrt(taux**2+tauy**2) + END DO + + RETURN +END SUBROUTINE ustarhb Index: LMDZ6/trunk/libf/phylmdiso/vdif_kcay.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/vdif_kcay.f90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/vdif_kcay.f90 (revision 5986) @@ -0,0 +1,674 @@ + +! $Header$ + +SUBROUTINE vdif_kcay(klon,klev,ngrid,dt, g, rconst, plev, temp, zlev, zlay, u, v, & + teta, cd, q2, q2diag, km, kn, ustar, l_mix) + + IMPLICIT NONE + + ! dt : pas de temps + ! g : g + ! zlev : altitude a chaque niveau (interface inferieure de la couche + ! de meme indice) + ! zlay : altitude au centre de chaque couche + ! u,v : vitesse au centre de chaque couche + ! (en entree : la valeur au debut du pas de temps) + ! teta : temperature potentielle au centre de chaque couche + ! (en entree : la valeur au debut du pas de temps) + ! cd : cdrag + ! (en entree : la valeur au debut du pas de temps) + ! q2 : $q^2$ au bas de chaque couche + ! (en entree : la valeur au debut du pas de temps) + ! (en sortie : la valeur a la fin du pas de temps) + ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque + ! couche) + ! (en sortie : la valeur a la fin du pas de temps) + ! kn : diffusivite turbulente des scalaires (au bas de chaque couche) + ! (en sortie : la valeur a la fin du pas de temps) + + ! ....................................................................... + INTEGER, INTENT(IN) :: klon,klev,ngrid + REAL,INTENT(IN) :: dt, g, rconst + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: plev + REAL,DIMENSION(klon,klev),INTENT(IN) :: temp + REAL,DIMENSION(klon),INTENT(IN) :: ustar + REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: zlev + REAL,DIMENSION(klon,klev),INTENT(IN) :: zlay + REAL,DIMENSION(klon,klev),INTENT(IN) :: u + REAL,DIMENSION(klon,klev),INTENT(IN) :: v + REAL,DIMENSION(klon,klev),INTENT(IN) :: teta + REAL,DIMENSION(klon),INTENT(IN) :: cd + REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2 + REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: q2diag + REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: km + REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: kn + INTEGER, INTENT(OUT) :: l_mix + + REAL,DIMENSION(klon) :: sq, sqz, long0 + REAL,DIMENSION(klon,klev+1) :: q2s,zz + REAL :: snstable,zq + + INTEGER :: iii + ! ....................................................................... + + ! nlay : nombre de couches + ! nlev : nombre de niveaux + ! ngrid : nombre de points de grille + ! unsdz : 1 sur l'epaisseur de couche + ! unsdzdec : 1 sur la distance entre le centre de la couche et le + ! centre de la couche inferieure + ! q : echelle de vitesse au bas de chaque couche + ! (valeur a la fin du pas de temps) + + ! ....................................................................... + INTEGER :: nlay, nlev + REAL, DIMENSION(klon,klev) :: unsdz + REAL, DIMENSION(klon, klev+1) :: unsdzdec,q + + ! ....................................................................... + + ! kmpre : km au debut du pas de temps + ! qcstat : q : solution stationnaire du probleme couple + ! (valeur a la fin du pas de temps) + ! q2cstat : q2 : solution stationnaire du probleme couple + ! (valeur a la fin du pas de temps) + + ! ....................................................................... + REAL, DIMENSION(klon, klev+1) :: kmpre + REAL :: qcstat + REAL :: q2cstat + REAL :: sss, sssq + ! ....................................................................... + + ! long : longueur de melange calculee selon Blackadar + + ! ....................................................................... + REAL, DIMENSION(klon, klev+1) :: long + ! ....................................................................... + + ! kmq3 : terme en q^3 dans le developpement de km + ! (valeur au debut du pas de temps) + ! kmcstat : valeur de km solution stationnaire du systeme {q2 ; du/dz} + ! (valeur a la fin du pas de temps) + ! knq3 : terme en q^3 dans le developpement de kn + ! mcstat : valeur de m solution stationnaire du systeme {q2 ; du/dz} + ! (valeur a la fin du pas de temps) + ! m2cstat : valeur de m2 solution stationnaire du systeme {q2 ; du/dz} + ! (valeur a la fin du pas de temps) + ! m : valeur a la fin du pas de temps + ! mpre : valeur au debut du pas de temps + ! m2 : valeur a la fin du pas de temps + ! n2 : valeur a la fin du pas de temps + + ! ....................................................................... + REAL :: kmq3 + REAL :: kmcstat + REAL :: knq3 + REAL :: mcstat + REAL :: m2cstat + REAL, DIMENSION(klon, klev+1) :: m,mpre,m2,n2 + ! ....................................................................... + + ! gn : intermediaire pour les coefficients de stabilite + ! gnmin : borne inferieure de gn (-0.23 ou -0.28) + ! gnmax : borne superieure de gn (0.0233) + ! gninf : vrai si gn est en dessous de sa borne inferieure + ! gnsup : vrai si gn est en dessus de sa borne superieure + ! gm : drole d'objet bien utile + ! ri : nombre de Richardson + ! sn : coefficient de stabilite pour n + ! snq2 : premier terme du developement limite de sn en q2 + ! sm : coefficient de stabilite pour m + ! smq2 : premier terme du developement limite de sm en q2 + + ! ....................................................................... + REAL :: gn,gm + REAL :: gnmin + REAL :: gnmax + LOGICAL :: gninf + LOGICAL :: gnsup + REAL, DIMENSION(klon, klev+1) :: sn, snq2, sm, smq2 + ! ....................................................................... + + ! kappa : consatnte de Von Karman (0.4) + ! long00 : longueur de reference pour le calcul de long (160) + ! a1,a2,b1,b2,c1 : constantes d'origine pour les coefficients + ! de stabilite (0.92/0.74/16.6/10.1/0.08) + ! cn1,cn2 : constantes pour sn + ! cm1,cm2,cm3,cm4 : constantes pour sm + + ! ....................................................................... + REAL :: kappa + REAL :: long00 + REAL :: a1, a2, b1, b2, c1 + REAL :: cn1, cn2 + REAL :: cm1, cm2, cm3, cm4 + ! ....................................................................... + + ! termq : termes en $q$ dans l'equation de q2 + ! termq3 : termes en $q^3$ dans l'equation de q2 + ! termqm2 : termes en $q*m^2$ dans l'equation de q2 + ! termq3m2 : termes en $q^3*m^2$ dans l'equation de q2 + + ! ....................................................................... + REAL :: termq + REAL :: termq3 + REAL :: termqm2 + REAL :: termq3m2 + ! ....................................................................... + + ! q2min : borne inferieure de q2 + ! q2max : borne superieure de q2 + + ! ....................................................................... + REAL :: q2min + REAL :: q2max + ! ....................................................................... + ! knmin : borne inferieure de kn + ! kmmin : borne inferieure de km + ! ....................................................................... + REAL :: knmin + REAL :: kmmin + ! ....................................................................... + INTEGER :: ilay, ilev, igrid + REAL :: tmp1, tmp2 + ! ....................................................................... + PARAMETER (kappa=0.4E+0) + PARAMETER (long00=160.E+0) + ! PARAMETER (gnmin=-10.E+0) + PARAMETER (gnmin=-0.28) + PARAMETER (gnmax=0.0233E+0) + PARAMETER (a1=0.92E+0) + PARAMETER (a2=0.74E+0) + PARAMETER (b1=16.6E+0) + PARAMETER (b2=10.1E+0) + PARAMETER (c1=0.08E+0) + PARAMETER (knmin=1.E-5) + PARAMETER (kmmin=1.E-5) + PARAMETER (q2min=1.E-5) + PARAMETER (q2max=1.E+2) + ! ym PARAMETER (nlay=klev) + ! ym PARAMETER (nlev=klev+1) + + PARAMETER (cn1=a2*(1.E+0-6.E+0*a1/b1)) + PARAMETER (cn2=-3.E+0*a2*(6.E+0*a1+b2)) + PARAMETER (cm1=a1*(1.E+0-3.E+0*c1-6.E+0*a1/b1)) + PARAMETER (cm2=a1*(-3.E+0*a2*((b2-3.E+0*a2)*(1.E+0-6.E+0*a1/b1)- & + 3.E+0*c1*(b2+6.E+0*a1)))) + PARAMETER (cm3=-3.E+0*a2*(6.E+0*a1+b2)) + PARAMETER (cm4=-9.E+0*a1*a2) + + LOGICAL :: first +! SAVE first +! DATA first/.TRUE./ +! !$OMP THREADPRIVATE(first) + ! ....................................................................... + ! traitment des valeur de q2 en entree + ! ....................................................................... + + ! Initialisation de q2 + nlay = klev + nlev = klev + 1 + +! Initialisation avec un schema d'equilibre +! CALL yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, cd, & +! q2diag, km, kn, ustar, l_mix) +! IF (first .AND. 1==1) THEN +! first = .FALSE. +! q2 = q2diag +! END IF +q2diag=0. + + DO ilev = 2, nlay + DO igrid = 1, ngrid + q2(igrid, ilev) = amax1(q2(igrid,ilev), q2min) + q(igrid, ilev) = sqrt(q2(igrid,ilev)) + END DO + END DO + + DO igrid = 1, ngrid + tmp1 = cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2) + q2(igrid, 1) = b1**(2.E+0/3.E+0)*tmp1 + q2(igrid, 1) = amax1(q2(igrid,1), q2min) + q(igrid, 1) = sqrt(q2(igrid,1)) + END DO + + ! ....................................................................... + ! les increments verticaux + ! ....................................................................... + + ! !!!!! allerte !!!!!c + ! !!!!! zlev n'est pas declare a nlev !!!!!c + ! !!!!! ----> + DO igrid = 1, ngrid + zlev(igrid, nlev) = zlay(igrid, nlay) + (zlay(igrid,nlay)-zlev(igrid,nlev & + -1)) + END DO + ! !!!!! <---- + ! !!!!! allerte !!!!!c + + DO ilay = 1, nlay + DO igrid = 1, ngrid + unsdz(igrid, ilay) = 1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay)) + END DO + END DO + DO igrid = 1, ngrid + unsdzdec(igrid, 1) = 1.E+0/(zlay(igrid,1)-zlev(igrid,1)) + END DO + DO ilay = 2, nlay + DO igrid = 1, ngrid + unsdzdec(igrid, ilay) = 1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1)) + END DO + END DO + DO igrid = 1, ngrid + unsdzdec(igrid, nlay+1) = 1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay)) + END DO + + ! ....................................................................... + ! le cisaillement et le gradient de temperature + ! ....................................................................... + + DO igrid = 1, ngrid + m2(igrid, 1) = (unsdzdec(igrid,1)*u(igrid,1))**2 + & + (unsdzdec(igrid,1)*v(igrid,1))**2 + m(igrid, 1) = sqrt(m2(igrid,1)) + mpre(igrid, 1) = m(igrid, 1) + END DO + + ! ----------------------------------------------------------------------- + DO ilev = 2, nlev - 1 + DO igrid = 1, ngrid + ! ----------------------------------------------------------------------- + + n2(igrid, ilev) = g*unsdzdec(igrid, ilev)*(teta(igrid,ilev)-teta(igrid, & + ilev-1))/(teta(igrid,ilev)+teta(igrid,ilev-1))*2.E+0 + ! n2(igrid,ilev)=0. + + ! ---> + ! on ne sais traiter que les cas stratifies. et l'ajustement + ! convectif est cense faire en sorte que seul des configurations + ! stratifiees soient rencontrees en entree de cette routine. + ! mais, bon ... on sait jamais (meme on sait que n2 prends + ! quelques valeurs negatives ... parfois) alors : + ! <--- + + IF (n2(igrid,ilev)<0.E+0) THEN + n2(igrid, ilev) = 0.E+0 + END IF + + m2(igrid, ilev) = (unsdzdec(igrid,ilev)*(u(igrid,ilev)-u(igrid, & + ilev-1)))**2 + (unsdzdec(igrid,ilev)*(v(igrid,ilev)-v(igrid, & + ilev-1)))**2 + m(igrid, ilev) = sqrt(m2(igrid,ilev)) + mpre(igrid, ilev) = m(igrid, ilev) + + ! ----------------------------------------------------------------------- + END DO + END DO + ! ----------------------------------------------------------------------- + + DO igrid = 1, ngrid + m2(igrid, nlev) = m2(igrid, nlev-1) + m(igrid, nlev) = m(igrid, nlev-1) + mpre(igrid, nlev) = m(igrid, nlev) + END DO + + ! ....................................................................... + ! calcul des fonctions de stabilite + ! ....................................................................... + + IF (l_mix==4) THEN + DO igrid = 1, ngrid + sqz(igrid) = 1.E-10 + sq(igrid) = 1.E-10 + END DO + DO ilev = 2, nlev - 1 + DO igrid = 1, ngrid + zq = sqrt(q2(igrid,ilev)) + sqz(igrid) = sqz(igrid) + zq*zlev(igrid, ilev)*(zlay(igrid,ilev)-zlay & + (igrid,ilev-1)) + sq(igrid) = sq(igrid) + zq*(zlay(igrid,ilev)-zlay(igrid,ilev-1)) + END DO + END DO + DO igrid = 1, ngrid + long0(igrid) = 0.2*sqz(igrid)/sq(igrid) + END DO + ELSE IF (l_mix==3) THEN + long0(igrid) = long00 + END IF + + ! (abd 5 2) print*,'LONG0=',long0 + + ! ----------------------------------------------------------------------- + DO ilev = 2, nlev - 1 + DO igrid = 1, ngrid + ! ----------------------------------------------------------------------- + + tmp1 = kappa*(zlev(igrid,ilev)-zlev(igrid,1)) + IF (l_mix>=10) THEN + long(igrid, ilev) = l_mix + ELSE + long(igrid, ilev) = tmp1/(1.E+0+tmp1/long0(igrid)) + END IF + long(igrid, ilev) = max(min(long(igrid,ilev),0.5*sqrt(q2(igrid,ilev))/ & + sqrt(max(n2(igrid,ilev),1.E-10))), 5.) + + gn = -long(igrid, ilev)**2/q2(igrid, ilev)*n2(igrid, ilev) + gm = long(igrid, ilev)**2/q2(igrid, ilev)*m2(igrid, ilev) + + gninf = .FALSE. + gnsup = .FALSE. + + IF (gngnmax) THEN + gnsup = .TRUE. + gn = gnmax + END IF + + sn(igrid, ilev) = cn1/(1.E+0+cn2*gn) + sm(igrid, ilev) = (cm1+cm2*gn)/((1.E+0+cm3*gn)*(1.E+0+cm4*gn)) + + IF ((gninf) .OR. (gnsup)) THEN + snq2(igrid, ilev) = 0.E+0 + smq2(igrid, ilev) = 0.E+0 + ELSE + snq2(igrid, ilev) = -gn*(-cn1*cn2/(1.E+0+cn2*gn)**2) + smq2(igrid, ilev) = -gn*(cm2*(1.E+0+cm3*gn)*(1.E+0+cm4*gn)-(cm3*( & + 1.E+0+cm4*gn)+cm4*(1.E+0+cm3*gn))*(cm1+cm2*gn))/((1.E+0+cm3*gn)*( & + 1.E+0+cm4*gn))**2 + END IF + + ! abd + ! if(ilev.le.57.and.ilev.ge.37) then + ! print*,'L=',ilev,' GN=',gn,' SM=',sm(igrid,ilev) + ! endif + ! ---> + ! la decomposition de Taylor en q2 n'a de sens que + ! dans les cas stratifies ou sn et sm sont quasi + ! proportionnels a q2. ailleurs on laisse le meme + ! algorithme car l'ajustement convectif fait le travail. + ! mais c'est delirant quand sn et snq2 n'ont pas le meme + ! signe : dans ces cas, on ne fait pas la decomposition. + ! <--- + + IF (snq2(igrid,ilev)*sn(igrid,ilev)<=0.E+0) snq2(igrid, ilev) = 0.E+0 + IF (smq2(igrid,ilev)*sm(igrid,ilev)<=0.E+0) smq2(igrid, ilev) = 0.E+0 + + ! Correction pour les couches stables. + ! Schema repris de JHoltzlag Boville, lui meme venant de... + + IF (1==1) THEN + snstable = 1. - zlev(igrid, ilev)/(700.*max(ustar(igrid),0.0001)) + snstable = 1. - zlev(igrid, ilev)/400. + snstable = max(snstable, 0.) + snstable = snstable*snstable + + ! abde print*,'SN ',ilev,sn(1,ilev),snstable + IF (sn(igrid,ilev) + DO ilev = 2, nlev - 1 + ! ----> + DO igrid = 1, ngrid + + ! ....................................................................... + + ! calcul des termes sources et puits de l'equation de q2 + ! ------------------------------------------------------ + + knq3 = kn(igrid, ilev)*snq2(igrid, ilev)/sn(igrid, ilev) + kmq3 = km(igrid, ilev)*smq2(igrid, ilev)/sm(igrid, ilev) + + termq = 0.E+0 + termq3 = 0.E+0 + termqm2 = 0.E+0 + termq3m2 = 0.E+0 + + tmp1 = dt*2.E+0*km(igrid, ilev)*m2(igrid, ilev) + tmp2 = dt*2.E+0*kmq3*m2(igrid, ilev) + termqm2 = termqm2 + dt*2.E+0*km(igrid, ilev)*m2(igrid, ilev) - & + dt*2.E+0*kmq3*m2(igrid, ilev) + termq3m2 = termq3m2 + dt*2.E+0*kmq3*m2(igrid, ilev) + + termq = termq - dt*2.E+0*kn(igrid, ilev)*n2(igrid, ilev) + & + dt*2.E+0*knq3*n2(igrid, ilev) + termq3 = termq3 - dt*2.E+0*knq3*n2(igrid, ilev) + + termq3 = termq3 - dt*2.E+0*q(igrid, ilev)**3/(b1*long(igrid,ilev)) + + ! ....................................................................... + + ! resolution stationnaire couplee avec le gradient de vitesse local + ! ----------------------------------------------------------------- + + ! -----{on cherche le cisaillement qui annule l'equation de q^2 + ! supposee en q3} + + tmp1 = termq + termq3 + tmp2 = termqm2 + termq3m2 + m2cstat = m2(igrid, ilev) - (tmp1+tmp2)/(dt*2.E+0*km(igrid,ilev)) + mcstat = sqrt(m2cstat) + + ! abde print*,'M2 L=',ilev,mpre(igrid,ilev),mcstat + + ! -----{puis on ecrit la valeur de q qui annule l'equation de m + ! supposee en q3} + + IF (ilev==2) THEN + kmcstat = 1.E+0/mcstat*(unsdz(igrid,ilev)*kmpre(igrid,ilev+1)*mpre( & + igrid,ilev+1)+unsdz(igrid,ilev-1)*cd(igrid)*(sqrt(u(igrid,3)**2+ & + v(igrid,3)**2)-mcstat/unsdzdec(igrid,ilev)-mpre(igrid, & + ilev+1)/unsdzdec(igrid,ilev+1))**2)/(unsdz(igrid,ilev)+unsdz(igrid, & + ilev-1)) + ELSE + kmcstat = 1.E+0/mcstat*(unsdz(igrid,ilev)*kmpre(igrid,ilev+1)*mpre( & + igrid,ilev+1)+unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)*mpre(igrid, & + ilev-1))/(unsdz(igrid,ilev)+unsdz(igrid,ilev-1)) + END IF + tmp2 = kmcstat/(sm(igrid,ilev)/q2(igrid,ilev))/long(igrid, ilev) + qcstat = tmp2**(1.E+0/3.E+0) + q2cstat = qcstat**2 + + ! ....................................................................... + + ! choix de la solution finale + ! --------------------------- + + q(igrid, ilev) = qcstat + q2(igrid, ilev) = q2cstat + m(igrid, ilev) = mcstat + ! abd if(ilev.le.57.and.ilev.ge.37) then + ! print*,'L=',ilev,' M2=',m2(igrid,ilev),m2cstat, + ! s 'N2=',n2(igrid,ilev) + ! abd endif + m2(igrid, ilev) = m2cstat + + ! ---> + ! pour des raisons simples q2 est minore + ! <--- + + IF (q2(igrid,ilev)gnmax) gn = gnmax + sn(igrid, ilev) = cn1/(1.E+0+cn2*gn) + sm(igrid, ilev) = (cm1+cm2*gn)/((1.E+0+cm3*gn)*(1.E+0+cm4*gn)) + kn(igrid, ilev) = long(igrid, ilev)*q(igrid, ilev)*sn(igrid, ilev) + km(igrid, ilev) = long(igrid, ilev)*q(igrid, ilev)*sm(igrid, ilev) + ! abd + ! if(ilev.le.57.and.ilev.ge.37) then + ! print*,'L=',ilev,' GN=',gn,' SM=',sm(igrid,ilev) + ! endif + + ! ....................................................................... + + END DO + + END DO + + ! ....................................................................... + + + DO igrid = 1, ngrid + kn(igrid, 1) = knmin + km(igrid, 1) = kmmin + ! kn(igrid,1)=cd(igrid) + ! km(igrid,1)=cd(igrid) + q2(igrid, nlev) = q2(igrid, nlev-1) + q(igrid, nlev) = q(igrid, nlev-1) + kn(igrid, nlev) = kn(igrid, nlev-1) + km(igrid, nlev) = km(igrid, nlev-1) + END DO + + ! CALCUL DE LA DIFFUSION VERTICALE DE Q2 + IF (1==1) THEN + + sss=0. + sssq=0. + ! WARNING : travail sur le point ig=1 ???? + DO ilev = 2, klev - 1 + sss = sss + plev(1, ilev-1) - plev(1, ilev+1) + sssq = sssq + (plev(1,ilev-1)-plev(1,ilev+1))*q2(1, ilev) + END DO + ! print*,'Q2moy avant',sssq/sss + ! print*,'Q2q20 ',(q2(1,ilev),ilev=1,10) + ! print*,'Q2km0 ',(km(1,ilev),ilev=1,10) + ! ! C'est quoi ca qu'etait dans l'original??? + ! do igrid=1,ngrid + ! q2(igrid,1)=10. + ! enddo + ! q2s=q2 + ! do iii=1,10 + ! call vdif_q2(dt,g,rconst,plev,temp,km,q2) + ! do ilev=1,klev+1 + ! write(iii+49,*) q2(1,ilev),zlev(1,ilev) + ! enddo + ! enddo + ! stop + ! do ilev=1,klev + ! print*,zlev(1,ilev),q2s(1,ilev),q2(1,ilev) + ! enddo + ! q2s=q2-q2s + ! do ilev=1,klev + ! print*,q2s(1,ilev),zlev(1,ilev) + ! enddo + DO ilev = 2, klev - 1 + sss = sss + plev(1, ilev-1) - plev(1, ilev+1) + sssq = sssq + (plev(1,ilev-1)-plev(1,ilev+1))*q2(1, ilev) + END DO + PRINT *, 'Q2moy apres', sssq/sss + + + DO ilev = 2, klev-1 + DO igrid = 1, ngrid + q2(igrid, ilev) = max(q2(igrid,ilev), q2min) + q(igrid, ilev) = sqrt(q2(igrid,ilev)) + + ! ....................................................................... + + ! calcul final de kn et km + ! ------------------------ + + gn = -long(igrid, ilev)**2/q2(igrid, ilev)*n2(igrid, ilev) + IF (gngnmax) gn = gnmax + sn(igrid, ilev) = cn1/(1.E+0+cn2*gn) + sm(igrid, ilev) = (cm1+cm2*gn)/((1.E+0+cm3*gn)*(1.E+0+cm4*gn)) + ! Correction pour les couches stables. + ! Schema repris de JHoltzlag Boville, lui meme venant de... + + IF (1==1) THEN + snstable = 1. - zlev(igrid, ilev)/(700.*max(ustar(igrid),0.0001)) + snstable = 1. - zlev(igrid, ilev)/400. + snstable = max(snstable, 0.) + snstable = snstable*snstable + + ! abde print*,'SN ',ilev,sn(1,ilev),snstable + IF (sn(igrid,ilev) the model starts with NP from start files created by ce0l + ! -> the model can run with longer time-steps. + ! 2016/11/30 (EV etienne.vignon@lmd.ipsl.fr) + ! On met tke (=q2/2) en entr??e plut??t que q2 + ! On corrige l'update de la tke + ! 2020/10/01 (EV) + ! On ajoute la dissipation de la TKE en diagnostique de sortie + ! + ! Inpout/Output : + !============== + ! tke : tke au bas de chaque couche + ! (en entree : la valeur au debut du pas de temps) + ! (en sortie : la valeur a la fin du pas de temps) + + ! Outputs: + !========== + ! eps: tke dissipation rate + ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque + ! couche) + ! (en sortie : la valeur a la fin du pas de temps) + ! kn : diffusivite turbulente des scalaires (au bas de chaque couche) + ! (en sortie : la valeur a la fin du pas de temps) + ! + !....................................................................... + + !======================================================================= + ! Declarations: + !======================================================================= + + + ! Inputs/Outputs + !---------------- + REAL dt, g, rconst + REAL plev(klon, klev+1), temp(klon, klev) + REAL ustar(klon) + REAL kmin, qmin, pblhmin(klon), coriol(klon) + REAL zlev(klon, klev+1) + REAL zlay(klon, klev) + REAL u(klon, klev) + REAL v(klon, klev) + REAL teta(klon, klev) + REAL cd(klon) + REAL tke(klon, klev+1) + REAL eps(klon,klev+1) + REAL unsdz(klon, klev) + REAL unsdzdec(klon, klev+1) + REAL kn(klon, klev+1) + REAL km(klon, klev+1) + INTEGER iflag_pbl, ngrid, nsrf + INTEGER ni(klon) + +!FC + REAL drgpro(klon,klev) + REAL winds(klon,klev) + + ! Local + !------- + + REAL q2(klon, klev+1) + REAL kmpre(klon, klev+1), tmp2, qpre + REAL mpre(klon, klev+1) + REAL kq(klon, klev+1) + REAL ff(klon, klev+1), delta(klon, klev+1) + REAL aa(klon, klev+1), aa0, aa1 + INTEGER nlay, nlev + + INTEGER ig, jg, k + REAL ri, zrif, zalpha, zsm, zsn + REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev) + REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1) + REAL dtetadz(klon, klev+1) + REAL m2cstat, mcstat, kmcstat + REAL l(klon, klev+1) + REAL zz(klon, klev+1) + INTEGER iter + REAL dissip(klon,klev), tkeprov,tkeexp, shear(klon,klev), buoy(klon,klev) + REAL :: disseff + + REAL,SAVE :: rifc + !$OMP THREADPRIVATE(rifc) + REAL,SAVE :: seuilsm, seuilalpha + !$OMP THREADPRIVATE(seuilsm, seuilalpha) + + REAL frif, falpha, fsm + REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), & + lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev) + + CHARACTER (len = 20) :: modname = 'yamada4' + CHARACTER (len = 80) :: abort_message + + + + ! Fonctions utilis??es + !-------------------- + + frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) + falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri) + fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri)) + + + + IF (new_yamada4) THEN +! Corrections et reglages issus du travail de these d'Etienne Vignon. + rifc=frif(ric) + seuilsm=fsm(frif(ric)) + seuilalpha=falpha(frif(ric)) + ELSE + rifc=0.191 + seuilalpha=1.12 + seuilsm=0.085 + ENDIF + +!=============================================================================== +! Flags, tests et ??valuations de constantes +!=============================================================================== + +! On utilise ou non la routine de Holstalg Boville pour les cas tres stables + + + IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN + abort_message='probleme de coherence dans appel a MY' + CALL abort_physic(modname,abort_message,1) + END IF + + + nlay = klev + nlev = klev + 1 + + +!======================================================================== +! Calcul des increments verticaux +!========================================================================= + + +! Attention: zlev n'est pas declare a nlev + DO ig = 1, ngrid + zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1)) + END DO + + + DO k = 1, nlay + DO ig = 1, ngrid + unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k)) + END DO + END DO + DO ig = 1, ngrid + unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1)) + END DO + DO k = 2, nlay + DO ig = 1, ngrid + unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1)) + END DO + END DO + DO ig = 1, ngrid + unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay)) + END DO + +!========================================================================= +! Richardson number and stability functions +!========================================================================= + +! initialize arrays: + + m2(1:ngrid, :) = 0.0 + sm(1:ngrid, :) = 0.0 + rif(1:ngrid, :) = 0.0 + +!------------------------------------------------------------ + DO k = 2, klev + + DO ig = 1, ngrid + dz(ig, k) = zlay(ig, k) - zlay(ig, k-1) + m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, & + k-1))**2)/(dz(ig,k)*dz(ig,k)) + dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k) + n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k)) + ri = n2(ig, k)/max(m2(ig,k), 1.E-10) + IF (ri0.) THEN + q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2 + ELSE + q2(ig, k) = (qpre/(1.-aa(ig,k)*qpre))**2 + END IF + ! else ! iflag_pbl=9 + ! if (aa(ig,k)*qpre.gt.0.9) then + ! q2(ig,k)=(qpre*10.)**2 + ! else + ! q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2 + ! endif + ! endif + q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4) + END DO + END DO + + ELSE IF (iflag_pbl>=10) THEN + + shear(:,:)=0. + buoy(:,:)=0. + dissip(:,:)=0. + km(:,:)=0. + + IF (yamada4_num>=1) THEN + + DO k = 2, klev - 1 + DO ig=1,ngrid + q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4) + km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k) + shear(ig,k)=km(ig, k)*m2(ig, k) + buoy(ig,k)=km(ig, k)*m2(ig, k)*(-1.*rif(ig,k)) +! dissip(ig,k)=min(max(((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)),1.E-12),1.E4) + dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)) + ENDDO + ENDDO + + IF (yamada4_num==1) THEN ! Schema du MAR tel quel + DO k = 2, klev - 1 + DO ig=1,ngrid + tkeprov=q2(ig,k)/ydeux + tkeprov= tkeprov* & + & (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ & + & (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k)+drgpro(ig,k)*tkeprov)) + q2(ig,k)=tkeprov*ydeux + ENDDO + ENDDO + ELSE IF (yamada4_num==2) THEN ! version modifiee avec integration exacte pour la dissipation + DO k = 2, klev - 1 + DO ig=1,ngrid + tkeprov=q2(ig,k)/ydeux + disseff=dissip(ig,k)-min(0.,buoy(ig,k)) + tkeprov = tkeprov/(1.+dt*disseff/(2.*tkeprov))**2 + tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))) + q2(ig,k)=tkeprov*ydeux + ! En cas stable, on traite la flotabilite comme la + ! dissipation, en supposant que buoy/q2^3 est constant. + ! Puis on prend la solution exacte + ENDDO + ENDDO + ELSE IF (yamada4_num==3) THEN ! version modifiee avec integration exacte pour la dissipation + DO k = 2, klev - 1 + DO ig=1,ngrid + tkeprov=q2(ig,k)/ydeux + disseff=dissip(ig,k)-min(0.,buoy(ig,k)) + tkeprov=tkeprov*exp(-dt*disseff/tkeprov) + tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))) + q2(ig,k)=tkeprov*ydeux + ! En cas stable, on traite la flotabilite comme la + ! dissipation, en supposant que buoy/q2^3 est constant. + ! Puis on prend la solution exacte + ENDDO + ENDDO + ELSE IF (yamada4_num==4) THEN ! version modifiee avec integration exacte pour la dissipation + DO k = 2, klev - 1 + DO ig=1,ngrid + tkeprov=q2(ig,k)/ydeux + tkeprov= tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))) + tkeprov= tkeprov* & + & tkeprov/ & + & (tkeprov+dt*((-1.)*min(0.,buoy(ig,k))+dissip(ig,k))) + q2(ig,k)=tkeprov*ydeux + ! En cas stable, on traite la flotabilite comme la + ! dissipation, en supposant que buoy/q2^3 est constant. + ! Puis on prend la solution exacte + ENDDO + ENDDO + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !! Attention, yamada4_num=5 est inexacte car néglige les termes de flottabilité + !! en conditions instables + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation + DO k = 2, klev - 1 + DO ig=1,ngrid + tkeprov=q2(ig,k)/ydeux +!FC on ajoute la dissipation due aux arbres + disseff=dissip(ig,k)-min(0.,buoy(ig,k)) + drgpro(ig,k)*tkeprov + tkeexp=exp(-dt*disseff/tkeprov) +! on prend en compte la tke cree par les arbres + winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) + tkeprov= (shear(ig,k)+ & + & drgpro(ig,k)*(winds(ig,k))**3)*tkeprov/disseff*(1.-tkeexp)+tkeprov*tkeexp + q2(ig,k)=tkeprov*ydeux + ! En cas stable, on traite la flotabilite comme la + ! dissipation, en supposant que buoy/q2^3 est constant. + ! Puis on prend la solution exacte + ENDDO + ENDDO + ELSE IF (yamada4_num==6) THEN ! version modifiee avec integration exacte pour la dissipation + DO k = 2, klev - 1 + DO ig=1,ngrid + ! En cas stable, on traite la flotabilite comme la + ! dissipation, en supposant que dissipeff/TKE est constant. + ! Puis on prend la solution exacte + ! With drag and dissipation from high vegetation (EV & FC, 05/10/2020) + winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) + tkeprov=q2(ig,k)/ydeux + tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3,0.)*dt + disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3) + drgpro(ig,k)*tkeprov + tkeexp=exp(-dt*disseff/tkeprov) + tkeprov= tkeprov*tkeexp + q2(ig,k)=tkeprov*ydeux + + ENDDO + ENDDO + ENDIF + + DO k = 2, klev - 1 + DO ig=1,ngrid + q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4) + ENDDO + ENDDO + + ELSE + + DO k = 2, klev - 1 + km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k) + q2(1:ngrid, k) = q2(1:ngrid, k) + ydeux*dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) +! q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k)) + q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4) + q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1)) +! q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1)) + q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k) + END DO + + ENDIF + + ELSE + abort_message='Cas nom prevu dans yamada4' + CALL abort_physic(modname,abort_message,1) + + END IF ! Fin du cas 8 + + + ! ==================================================================== + ! Calcul des coefficients de melange + ! ==================================================================== + + DO k = 2, klev + DO ig = 1, ngrid + zq = sqrt(q2(ig,k)) + km(ig, k) = l(ig, k)*zq*sm(ig, k) ! For momentum + kn(ig, k) = km(ig, k)*alpha(ig, k) ! For scalars + kq(ig, k) = l(ig, k)*zq*0.2 ! For TKE + END DO + END DO + + + !==================================================================== + ! Transport diffusif vertical de la TKE par la TKE + !==================================================================== + + + ! initialize near-surface and top-layer mixing coefficients + !........................................................... + + kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface + kq(1:ngrid, klev+1) = 0 ! zero at the top + + ! Transport diffusif vertical de la TKE. + !....................................... + + IF (iflag_vdif_q2==1) THEN + q2(1:ngrid, 1) = q2(1:ngrid, 2) + CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2) + END IF + + + !==================================================================== + ! Traitement particulier pour les cas tres stables, introduction d'une + ! longueur de m??lange minimale + !==================================================================== + ! + ! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model + ! Holtslag A.A.M. and Boville B.A. + ! J. Clim., 6, 1825-1842, 1993 + + + IF (hboville) THEN + + + IF (prt_level>1) THEN + WRITE(lunout,*) 'YAMADA4 0' + END IF + + DO ig = 1, ngrid + coriol(ig) = 1.E-4 + pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5) + END DO + + IF (1==1) THEN + IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN + + DO k = 2, klev + DO ig = 1, ngrid + IF (teta(ig,2)>teta(ig,1)) THEN + qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2 + kmin = kap*zlev(ig, k)*qmin + ELSE + kmin = -1. ! kmin n'est utilise que pour les SL stables. + END IF + IF (kn(ig,k)teta(ig,1)) THEN + qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2 + kmin = kap*zlev(ig, k)*qmin + ELSE + kmin = -1. ! kmin n'est utilise que pour les SL stables. + END IF + IF (kn(ig,k)1) THEN + WRITE(lunout,*)'YAMADA4 1' + END IF !(prt_level>1) THEN + + + !====================================================== + ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane + !====================================================== + ! + ! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer + ! Abdella K and McFarlane N + ! J. Atmos. Sci., 54, 1850-1867, 1997 + + ! Diagnostique pour stokage + !.......................... + + IF (1==0) THEN + rino = rif + smyam(1:ngrid, 1) = 0. + styam(1:ngrid, 1) = 0. + lyam(1:ngrid, 1) = 0. + knyam(1:ngrid, 1) = 0. + w2yam(1:ngrid, 1) = 0. + t2yam(1:ngrid, 1) = 0. + + smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev) + styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev) + lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev) + knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev) + + + ! Calcul de w'2 et T'2 + !....................... + + w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + & + lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ & + sqrt(q2(1:ngrid,2:klev)) + + t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* & + dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* & + lyam(1:ngrid, 2:klev) + END IF + + + +!============================================================================ +! Mise a jour de la tke +!============================================================================ + + IF (new_yamada4) THEN + DO k=1,klev+1 + tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux + ENDDO + ELSE + DO k=1,klev+1 + tke(1:ngrid,k)=q2(1:ngrid,k) + ENDDO + ENDIF + + +!============================================================================ +! Diagnostique de la dissipation et vitesse verticale +!============================================================================ + +! Diagnostics + + eps(:,:)=0. + wprime(1:ngrid,:,nsrf)=0. + DO k=2,klev + DO ig=1,ngrid + eps(ig,k)=dissip(ig,k) + jg=ni(ig) + wprime(jg,k,nsrf)=sqrt(MAX(1./3*q2(ig,k),0.)) + ENDDO + ENDDO + + +!============================================================================= + + RETURN + + +END SUBROUTINE yamada4 + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2) + + USE dimphy, only : klev,klon + IMPLICIT NONE + +! vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE +! avec un schema implicite en temps avec +! inversion d'un syst??me tridiagonal +! +! Reference: Description of the interface with the surface and +! the computation of the turbulet diffusion in LMDZ +! Technical note on LMDZ +! Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y +! +!============================================================================ +! Declarations +!============================================================================ + + REAL plev(klon, klev+1) + REAL temp(klon, klev) + REAL timestep + REAL gravity, rconst + REAL kstar(klon, klev+1), zz + REAL kmy(klon, klev+1) + REAL q2(klon, klev+1) + REAL deltap(klon, klev+1) + REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1) + INTEGER ngrid + + INTEGER i, k + + +!========================================================================= +! Calcul +!========================================================================= + + DO k = 1, klev + DO i = 1, ngrid + zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k)) + kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ & + (plev(i,k)-plev(i,k+1))*timestep + END DO + END DO + + DO k = 2, klev + DO i = 1, ngrid + deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1)) + END DO + END DO + DO i = 1, ngrid + deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2)) + deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1)) + denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev) + alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1) + beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1) + END DO + + DO k = klev, 2, -1 + DO i = 1, ngrid + denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + & + kstar(i, k-1) + alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k) + beta(i, k) = kstar(i, k-1)/denom(i, k) + END DO + END DO + + ! Si on recalcule q2(1) + !....................... + IF (1==0) THEN + DO i = 1, ngrid + denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1) + q2(i, 1) = (q2(i,1)*deltap(i,1)+kstar(i,1)*alpha(i,2))/denom(i, 1) + END DO + END IF + + + DO k = 2, klev + 1 + DO i = 1, ngrid + q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1) + END DO + END DO + + RETURN +END SUBROUTINE vdif_q2 +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2) + + USE dimphy, only : klev,klon + IMPLICIT NONE + +! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE +! avec un schema explicite en temps + + +!==================================================== +! Declarations +!==================================================== + + REAL plev(klon, klev+1) + REAL temp(klon, klev) + REAL timestep + REAL gravity, rconst + REAL kstar(klon, klev+1), zz + REAL kmy(klon, klev+1) + REAL q2(klon, klev+1) + REAL deltap(klon, klev+1) + REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1) + INTEGER ngrid + INTEGER i, k + + +!================================================== +! Calcul +!================================================== + + DO k = 1, klev + DO i = 1, ngrid + zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k)) + kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ & + (plev(i,k)-plev(i,k+1))*timestep + END DO + END DO + + DO k = 2, klev + DO i = 1, ngrid + deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1)) + END DO + END DO + DO i = 1, ngrid + deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2)) + deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1)) + END DO + + DO k = klev, 2, -1 + DO i = 1, ngrid + q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, & + k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k) + END DO + END DO + + DO i = 1, ngrid + q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1) + q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, & + klev)))/deltap(i, klev+1) + END DO + + RETURN +END SUBROUTINE vdif_q2e + +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +SUBROUTINE mixinglength(ni, nsrf, ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, lmix) + + + + USE dimphy, only : klev,klon + USE yamada_ini_mod, only : l0 + USE phys_state_var_mod, only: zstd, zsig, zmea + USE phys_local_var_mod, only: l_mixmin, l_mix + USE yamada_ini_mod, only : kap, kapb + + ! zstd: ecart type de la'altitud e sous-maille + ! zmea: altitude moyenne sous maille + ! zsig: pente moyenne de le maille + + USE geometry_mod, only: cell_area + ! aire_cell: aire de la maille + + IMPLICIT NONE +!************************************************************************* +! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence +! avec la formule de Blackadar +! Calcul d'un minimum en fonction de l'orographie sous-maille: +! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange +! induit par les circulations meso et submeso ??chelles. +! +! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere +! Blackadar A.K. +! J. Geophys. Res., 64, No 8, 1962 +! +! * An evaluation of neutral and convective planetary boundary-layer parametrisations relative +! to large eddy simulations +! Ayotte K et al +! Boundary Layer Meteorology, 79, 131-175, 1996 +! +! +! * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts +! Van de Wiel B.J.H et al +! Boundary-Lay Meteorol, 128, 103-166, 2008 +! +! +! Histoire: +!---------- +! * premi??re r??daction, Etienne et Frederic, 09/06/2016 +! +! *********************************************************************** + +!================================================================== +! Declarations +!================================================================== + +! Inputs +!------- + INTEGER ni(klon) ! indice sur la grille original (non restreinte) + INTEGER nsrf ! Type de surface + INTEGER ngrid ! Nombre de points concern??s sur l'horizontal + INTEGER iflag_pbl ! Choix du sch??ma de turbulence + REAL pbl_lmixmin_alpha ! on active ou non le calcul de la longueur de melange minimum en fonction du relief + REAL lmixmin ! Minimum absolu de la longueur de m??lange + REAL zlay(klon, klev) ! altitude du centre de la couche + REAL zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche + REAL u(klon, klev) ! vitesse du vent zonal + REAL v(klon, klev) ! vitesse du vent meridional + REAL q2(klon, klev+1) ! energie cin??tique turbulente + REAL n2(klon, klev+1) ! frequence de Brunt-Vaisala + +!In/out +!------- + +! Outputs +!--------- + + REAL lmix(klon, klev+1) ! Longueur de melange + + +! Local +!------- + + INTEGER ig,jg, k + REAL h_oro(klon) + REAL hlim(klon) + REAL zq + REAL sq(klon), sqz(klon) + REAL fl, zzz, zl0, zq2, zn2 + REAL famorti, zzzz, zh_oro, zhlim + REAL l1(klon, klev+1), l2(klon,klev+1) + REAL winds(klon, klev) + REAL xcell + REAL zstdslope(klon) + REAL lmax + REAL l2strat, l2neutre, extent + REAL l2limit(klon) +!=============================================================== +! Fonctions utiles +!=============================================================== + +! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996 +!.......................................................................... + + fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, & + k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( & + max(n2(ig,k),1.E-10))), 1.E-5) + +! Fonction d'amortissement de la turbulence au dessus de la montagne +! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code +!..................................................................... + + famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1. + + IF (ngrid==0) RETURN + + +!===================================================================== +! CALCUL de la LONGUEUR de m??lange suivant BLACKADAR: l1 +!===================================================================== + + l1(1:ngrid,:)=0. + IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN + + + ! Iterative computation of l0 + ! This version is kept for iflag_pbl only for convergence + ! with NPv3.1 Cmip5 simulations + !................................................................... + + DO ig = 1, ngrid + sq(ig) = 1.E-10 + sqz(ig) = 1.E-10 + END DO + DO k = 2, klev - 1 + DO ig = 1, ngrid + zq = sqrt(q2(ig,k)) + sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1)) + sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1)) + END DO + END DO + DO ig = 1, ngrid + l0(ig) = 0.2*sqz(ig)/sq(ig) + END DO + DO k = 2, klev + DO ig = 1, ngrid + l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k)) + END DO + END DO + + ELSE + + + ! In all other case, the assymptotic mixing length l0 is imposed (150m) + !...................................................................... + + l0(1:ngrid) = 150. + DO k = 2, klev + DO ig = 1, ngrid + l1(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k)) + END DO + END DO + + END IF + +!=========================================================================================== +! CALCUL d'une longueur de melange minimum en fonctions de la topographie sous maille: l2 +! si pbl_lmixmin_alpha=TRUE et si on se trouve sur de la terre ( pas actif sur les +! glacier, la glace de mer et les oc??ans) +!=========================================================================================== + + l2(1:ngrid,:)=0.0 + l_mixmin(1:ngrid,:,nsrf)=0. + l_mix(1:ngrid,:,nsrf)=0. + hlim(1:ngrid)=0. + + IF (nsrf .EQ. 1) THEN + +! coefficients +!-------------- + + extent=2. ! On ??tend l'impact du relief jusqu'?? extent*h, extent >1. + lmax=150. ! Longueur de m??lange max dans l'absolu + +! calculs +!--------- + + DO ig=1,ngrid + + ! On calcule la hauteur du relief + !................................. + ! On ne peut pas prendre zstd seulement pour caracteriser le relief sous maille + ! car sur un terrain pentu mais sans relief, zstd est non nul (comme en Antarctique, C. Genthon) + ! On corrige donc zstd avec l'ecart type de l'altitude dans la maille sans relief + ! (en gros, une maille de taille xcell avec une pente constante zstdslope) + jg=ni(ig) +! IF (zsig(jg) .EQ. 0.) THEN +! zstdslope(ig)=0. +! ELSE +! xcell=sqrt(cell_area(jg)) +! zstdslope(ig)=max((xcell*zsig(jg)-zmea(jg))**3 /(3.*zsig(jg)),0.) +! zstdslope(ig)=sqrt(zstdslope(ig)) +! END IF + +! h_oro(ig)=max(zstd(jg)-zstdslope(ig),0.) ! Hauteur du relief + h_oro(ig)=zstd(jg) + hlim(ig)=extent*h_oro(ig) + ENDDO + + l2limit(1:ngrid)=0. + + DO k=2,klev + DO ig=1,ngrid + winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) + IF (zlev(ig,k) .LE. h_oro(ig)) THEN ! sous l'orographie + l2strat= kapb*pbl_lmixmin_alpha*winds(ig,k)/sqrt(max(n2(ig,k),1.E-10)) ! si stratifi??, amplitude d'oscillation * kappab (voir Van de Wiel et al 2008) + l2neutre=kap*zlev(ig,k)*h_oro(ig)/(kap*zlev(ig,k)+h_oro(ig)) ! Dans le cas neutre, formule de blackadar. tend asymptotiquement vers h + l2neutre=MIN(l2neutre,lmax) ! On majore par lmax + l2limit(ig)=MIN(l2neutre,l2strat) ! Calcule de l2 (minimum de la longueur en cas neutre et celle en situation stratifi??e) + l2(ig,k)=l2limit(ig) + + ELSE IF (zlev(ig,k) .LE. hlim(ig)) THEN ! Si on est au dessus des montagnes, mais affect?? encore par elles + + ! Au dessus des montagnes, on prend la l2limit au sommet des montagnes + ! (la derni??re calcul??e dans la boucle k, vu que k est un indice croissant avec z) + ! et on multiplie l2limit par une fonction qui d??croit entre h et hlim + l2(ig,k)=l2limit(ig)*famorti(zlev(ig,k),h_oro(ig), hlim(ig)) + ELSE ! Au dessus de extent*h, on prend l2=l0 + l2(ig,k)=0. + END IF + ENDDO + ENDDO + ENDIF ! pbl_lmixmin_alpha + +!================================================================================== +! On prend le max entre la longueur de melange de blackadar et celle calcul??e +! en fonction de la topographie +!=================================================================================== + + + DO k=1,klev+1 + DO ig=1,ngrid + lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin) + ENDDO + ENDDO + +! Diagnostics + + DO k=1,klev+1 + DO ig=1,ngrid + jg=ni(ig) + l_mix(jg,k,nsrf)=lmix(ig,k) + l_mixmin(jg,k,nsrf)=MAX(l2(ig,k),lmixmin) + ENDDO + ENDDO + + + +END SUBROUTINE mixinglength Index: LMDZ6/trunk/libf/phylmdiso/yamada_c.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/yamada_c.F90 (revision 5986) +++ LMDZ6/trunk/libf/phylmdiso/yamada_c.F90 (revision 5986) @@ -0,0 +1,481 @@ +! +! $Header$ +! + SUBROUTINE yamada_c(ngrid,timestep,plev,play & + & ,pu,pv,pt,d_u,d_v,d_t,cd,q2,km,kn,kq,d_t_diss,ustar & + & ,iflag_pbl) + USE dimphy, ONLY: klon, klev + USE print_control_mod, ONLY: prt_level + USE ioipsl_getin_p_mod, ONLY : getin_p + + USE yomcst_mod_h +IMPLICIT NONE + +! +! timestep : pas de temps +! g : g +! zlev : altitude a chaque niveau (interface inferieure de la couche +! de meme indice) +! zlay : altitude au centre de chaque couche +! u,v : vitesse au centre de chaque couche +! (en entree : la valeur au debut du pas de temps) +! teta : temperature potentielle au centre de chaque couche +! (en entree : la valeur au debut du pas de temps) +! cd : cdrag +! (en entree : la valeur au debut du pas de temps) +! q2 : $q^2$ au bas de chaque couche +! (en entree : la valeur au debut du pas de temps) +! (en sortie : la valeur a la fin du pas de temps) +! km : diffusivite turbulente de quantite de mouvement (au bas de chaque +! couche) +! (en sortie : la valeur a la fin du pas de temps) +! kn : diffusivite turbulente des scalaires (au bas de chaque couche) +! (en sortie : la valeur a la fin du pas de temps) +! +! iflag_pbl doit valoir entre 6 et 9 +! l=6, on prend systematiquement une longueur d'equilibre +! iflag_pbl=6 : MY 2.0 +! iflag_pbl=7 : MY 2.0.Fournier +! iflag_pbl=8/9 : MY 2.5 +! iflag_pbl=8 with special obsolete treatments for convergence +! with Cmpi5 NPv3.1 simulations +! iflag_pbl=10/11 : New scheme M2 and N2 explicit and dissiptation exact +! iflag_pbl=12 = 11 with vertical diffusion off q2 +! +! 2013/04/01 (FH hourdin@lmd.jussieu.fr) +! Correction for very stable PBLs (iflag_pbl=10 and 11) +! iflag_pbl=8 converges numerically with NPv3.1 +! iflag_pbl=11 -> the model starts with NP from start files created by ce0l +! -> the model can run with longer time-steps. +!....................................................................... + + REAL, DIMENSION(klon,klev) :: d_u,d_v,d_t + REAL, DIMENSION(klon,klev) :: pu,pv,pt + REAL, DIMENSION(klon,klev) :: d_t_diss + + REAL timestep + real plev(klon,klev+1) + real play(klon,klev) + real ustar(klon) + real kmin,qmin,pblhmin(klon),coriol(klon) + REAL zlev(klon,klev+1) + REAL zlay(klon,klev) + REAL zu(klon,klev) + REAL zv(klon,klev) + REAL zt(klon,klev) + REAL teta(klon,klev) + REAL cd(klon) + REAL q2(klon,klev+1),qpre + REAL unsdz(klon,klev) + REAL unsdzdec(klon,klev+1) + + REAL km(klon,klev) + REAL kmpre(klon,klev+1),tmp2 + REAL mpre(klon,klev+1) + REAL kn(klon,klev) + REAL kq(klon,klev) + real ff(klon,klev+1),delta(klon,klev+1) + real aa(klon,klev+1),aa0,aa1 + integer iflag_pbl,ngrid + integer nlay,nlev + + logical first + integer ipas + save first,ipas +!FH/IM data first,ipas/.true.,0/ + data first,ipas/.false.,0/ +!$OMP THREADPRIVATE( first,ipas) + INTEGER, SAVE :: iflag_tke_diff=0 +!$OMP THREADPRIVATE(iflag_tke_diff) + + + integer ig,k + + + real ri,zrif,zalpha,zsm,zsn + real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev) + + real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1) + REAL, DIMENSION(klon,klev+1) :: km2,kn2,sqrtq + real dtetadz(klon,klev+1) + real m2cstat,mcstat,kmcstat + real l(klon,klev+1) + real leff(klon,klev+1) + real,allocatable,save :: l0(:) +!$OMP THREADPRIVATE(l0) + real sq(klon),sqz(klon),zz(klon,klev+1) + integer iter + + real ric,rifc,b1,kap + save ric,rifc,b1,kap + data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/ +!$OMP THREADPRIVATE(ric,rifc,b1,kap) + real frif,falpha,fsm + real fl,zzz,zl0,zq2,zn2 + + real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev) + real lyam(klon,klev),knyam(klon,klev) + real w2yam(klon,klev),t2yam(klon,klev) + logical,save :: firstcall=.true. +!$OMP THREADPRIVATE(firstcall) + CHARACTER(len=20),PARAMETER :: modname="yamada_c" +REAL, DIMENSION(klon,klev+1) :: fluxu,fluxv,fluxt +REAL, DIMENSION(klon,klev+1) :: dddu,dddv,dddt +REAL, DIMENSION(klon,klev) :: exner,masse +REAL, DIMENSION(klon,klev+1) :: masseb,q2old,q2neg + LOGICAL okiophys + + frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) + falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri) + fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri)) + fl(zzz,zl0,zq2,zn2)= & + & max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig)) & + & ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.) + + + okiophys=klon==1 + if (firstcall) then + CALL getin_p('iflag_tke_diff',iflag_tke_diff) + allocate(l0(klon)) + firstcall=.false. + endif + + IF (ngrid<=0) RETURN ! Bizarre : on n a pas ce probeleme pour coef_diff_turb + + nlay=klev + nlev=klev+1 + + +!------------------------------------------------------------------------- +! Computation of conservative source terms from the turbulent tendencies +!------------------------------------------------------------------------- + + + zalpha=0.5 ! Anciennement 0.5. Essayer de voir pourquoi ? + zu(:,:)=pu(:,:)+zalpha*d_u(:,:) + zv(:,:)=pv(:,:)+zalpha*d_v(:,:) + zt(:,:)=pt(:,:)+zalpha*d_t(:,:) + + do k=1,klev + exner(:,k)=(play(:,k)/plev(:,1))**RKAPPA + masse(:,k)=(plev(:,k)-plev(:,k+1))/RG + teta(:,k)=zt(:,k)/exner(:,k) + enddo + +! Atmospheric mass at layer interfaces, where the TKE is computed + masseb(:,:)=0. + do k=1,klev + masseb(:,k)=masseb(:,k)+masse(:,k) + masseb(:,k+1)=masseb(:,k+1)+masse(:,k) + enddo + masseb(:,:)=0.5*masseb(:,:) + + zlev(:,1)=0. + zlay(:,1)=RCPD*teta(:,1)*(1.-exner(:,1)) + do k=1,klev-1 + zlay(:,k+1)=zlay(:,k)+0.5*RCPD*(teta(:,k)+teta(:,k+1))*(exner(:,k)-exner(:,k+1))/RG + zlev(:,k)=0.5*(zlay(:,k)+zlay(:,k+1)) ! PASBO + enddo + + fluxu(:,klev+1)=0. + fluxv(:,klev+1)=0. + fluxt(:,klev+1)=0. + + do k=klev,1,-1 + fluxu(:,k)=fluxu(:,k+1)+masse(:,k)*d_u(:,k) + fluxv(:,k)=fluxv(:,k+1)+masse(:,k)*d_v(:,k) + fluxt(:,k)=fluxt(:,k+1)+masse(:,k)*d_t(:,k)/exner(:,k) ! Flux de theta + enddo + + dddu(:,1)=2*zu(:,1)*fluxu(:,1) + dddv(:,1)=2*zv(:,1)*fluxv(:,1) + dddt(:,1)=(exner(:,1)-1.)*fluxt(:,1) + + do k=2,klev + dddu(:,k)=(zu(:,k)-zu(:,k-1))*fluxu(:,k) + dddv(:,k)=(zv(:,k)-zv(:,k-1))*fluxv(:,k) + dddt(:,k)=(exner(:,k)-exner(:,k-1))*fluxt(:,k) + enddo + dddu(:,klev+1)=0. + dddv(:,klev+1)=0. + dddt(:,klev+1)=0. + +#ifdef IOPHYS +if (okiophys) then + call iophys_ecrit('zlay',klev,'Geop','m',zlay) + call iophys_ecrit('teta',klev,'teta','K',teta) + call iophys_ecrit('temp',klev,'temp','K',zt) + call iophys_ecrit('pt',klev,'temp','K',pt) + call iophys_ecrit('pu',klev,'u','m/s',pu) + call iophys_ecrit('pv',klev,'v','m/s',pv) + call iophys_ecrit('d_u',klev,'d_u','m/s2',d_u) + call iophys_ecrit('d_v',klev,'d_v','m/s2',d_v) + call iophys_ecrit('d_t',klev,'d_t','K/s',d_t) + call iophys_ecrit('exner',klev,'exner','',exner) + call iophys_ecrit('masse',klev,'masse','',masse) + call iophys_ecrit('masseb',klev,'masseb','',masseb) +endif +#endif + + + + ipas=ipas+1 + + +!....................................................................... +! les increments verticaux +!....................................................................... +! +!!!!!! allerte !!!!!c +!!!!!! zlev n'est pas declare a nlev !!!!!c +!!!!!! ----> + DO ig=1,ngrid + zlev(ig,nlev)=zlay(ig,nlay) & + & +( zlay(ig,nlay) - zlev(ig,nlev-1) ) + ENDDO +!!!!!! <---- +!!!!!! allerte !!!!!c +! + DO k=1,nlay + DO ig=1,ngrid + unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k)) + ENDDO + ENDDO + DO ig=1,ngrid + unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1)) + ENDDO + DO k=2,nlay + DO ig=1,ngrid + unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1)) + ENDDO + ENDDO + DO ig=1,ngrid + unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay)) + ENDDO +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Computing M^2, N^2, Richardson numbers, stability functions +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + do k=2,klev + do ig=1,ngrid + dz(ig,k)=zlay(ig,k)-zlay(ig,k-1) + m2(ig,k)=((zu(ig,k)-zu(ig,k-1))**2+(zv(ig,k)-zv(ig,k-1))**2)/(dz(ig,k)*dz(ig,k)) + dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k) + n2(ig,k)=RG*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k)) +! n2(ig,k)=0. + ri=n2(ig,k)/max(m2(ig,k),1.e-10) + if (ri.lt.ric) then + rif(ig,k)=frif(ri) + else + rif(ig,k)=rifc + endif + if(rif(ig,k)<0.16) then + alpha(ig,k)=falpha(rif(ig,k)) + sm(ig,k)=fsm(rif(ig,k)) + else + alpha(ig,k)=1.12 + sm(ig,k)=0.085 + endif + zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k) + enddo + enddo + + + +!==================================================================== +! Computing the mixing length +!==================================================================== + +! Mise a jour de l0 + if (iflag_pbl==8.or.iflag_pbl==10) then + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Iterative computation of l0 +! This version is kept for iflag_pbl only for convergence +! with NPv3.1 Cmip5 simulations +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + do ig=1,ngrid + sq(ig)=1.e-10 + sqz(ig)=1.e-10 + enddo + do k=2,klev-1 + do ig=1,ngrid + zq=sqrt(q2(ig,k)) + sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1)) + sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1)) + enddo + enddo + do ig=1,ngrid + l0(ig)=0.2*sqz(ig)/sq(ig) + enddo + do k=2,klev + do ig=1,ngrid + l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k)) + enddo + enddo +! print*,'L0 cas 8 ou 10 ',l0 + + else + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! In all other case, the assymptotic mixing length l0 is imposed (100m) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + l0(:)=150. + do k=2,klev + do ig=1,ngrid + l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k)) + enddo + enddo +! print*,'L0 cas autres ',l0 + + endif + + +#ifdef IOPHYS +if (okiophys) then +call iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev)) +call iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev)) +call iophys_ecrit('Km2app',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev)) +call iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev)) +endif +#endif + + +IF (iflag_pbl<20) then + ! For diagnostics only + RETURN + +ELSE + +! print*,'OK1' + +! Evolution of TKE under source terms K M2 and K N2 + leff(:,:)=max(l(:,:),1.) + +!################################################################## +!# IF (iflag_pbl==29) THEN +!# STOP'Ne pas utiliser iflag_pbl=29' +!# km2(:,:)=km(:,:)*m2(:,:) +!# kn2(:,:)=kn2(:,:)*rif(:,:) +!# ELSEIF (iflag_pbl==25) THEN +! VERSION AVEC LA TKE EN MILIEU DE COUCHE +!# STOP'Ne pas utiliser iflag_pbl=25' +!# DO k=1,klev +!# km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) & +!# & /(masse(:,k)*timestep) +!# kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep) +!# leff(:,k)=0.5*(leff(:,k)+leff(:,k+1)) +!# ENDDO +!# km2(:,klev+1)=0. ; kn2(:,klev+1)=0. +!# ELSE +!################################################################# + + km2(:,:)=-(dddu(:,:)+dddv(:,:))/(masseb(:,:)*timestep) + kn2(:,:)=rcpd*dddt(:,:)/(masseb(:,:)*timestep) +! ENDIF + q2neg(:,:)=q2(:,:)+timestep*(km2(:,:)-kn2(:,:)) + q2(:,:)=min(max(q2neg(:,:),1.e-10),1.e4) + + +#ifdef IOPHYS +if (okiophys) then + call iophys_ecrit('km2',klev,'m2 conserv','m/s',km2(:,1:klev)) + call iophys_ecrit('kn2',klev,'n2 conserv','m/s',kn2(:,1:klev)) +endif +#endif + +! Dissipation of TKE + q2old(:,:)=q2(:,:) + q2(:,:)=1./(1./sqrt(q2(:,:))+timestep/(2*leff(:,:)*b1)) + q2(:,:)=q2(:,:)*q2(:,:) +! IF (iflag_pbl<=24) THEN + DO k=1,klev + d_t_diss(:,k)=(masseb(:,k)*(q2neg(:,k)-q2(:,k))+masseb(:,k+1)*(q2neg(:,k+1)-q2(:,k+1)))/(2.*rcpd*masse(:,k)) + ENDDO + +!################################################################### +! ELSE IF (iflag_pbl<=27) THEN +! DO k=1,klev +! d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd +! ENDDO +! ENDIF +! print*,'iflag_pbl ',d_t_diss +!################################################################### + + +! Compuation of stability functions +! IF (iflag_pbl/=29) THEN + DO k=1,klev + DO ig=1,ngrid + IF (ABS(km2(ig,k))<=1.e-20) THEN + rif(ig,k)=0. + ELSE + rif(ig,k)=min(kn2(ig,k)/km2(ig,k),rifc) + ENDIF + IF (rif(ig,k).lt.0.16) THEN + alpha(ig,k)=falpha(rif(ig,k)) + sm(ig,k)=fsm(rif(ig,k)) + else + alpha(ig,k)=1.12 + sm(ig,k)=0.085 + endif + ENDDO + ENDDO +! ENDIF + +! Computation of turbulent diffusivities +! IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN +! DO k=2,klev +! sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1))) +! ENDDO +! ELSE + kq(:,:)=0. + DO k=1,klev + ! Coefficient au milieu des couches pour diffuser la TKE + kq(:,k)=0.5*leff(:,k)*sqrt(q2(:,k))*0.2 + ENDDO + +#ifdef IOPHYS +if (okiophys) then +call iophys_ecrit('q2b',klev,'KTE inter','m2/s',q2(:,1:klev)) +endif +#endif + + IF (iflag_tke_diff==1) THEN + CALL vdif_q2(timestep, RG, RD, ngrid, plev, pt, kq, q2) + ENDIF + + km(:,:)=0. + kn(:,:)=0. + DO k=1,klev + km(:,k)=leff(:,k)*sqrt(q2(:,k))*sm(:,k) + kn(:,k)=km(:,k)*alpha(:,k) + ENDDO + + +#ifdef IOPHYS +if (okiophys) then +call iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev)) +call iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev)) +call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev)) +call iophys_ecrit('q2neg',klev,'KTE non bornee','m2/s',q2neg(:,1:klev)) +call iophys_ecrit('alpha',klev,'alpha','',alpha(:,1:klev)) +call iophys_ecrit('sm',klev,'sm','',sm(:,1:klev)) +call iophys_ecrit('q2f',klev,'KTE finale','m2/s',q2(:,1:klev)) +call iophys_ecrit('kmf',klev,'Kz final','m2/s',km(:,1:klev)) +call iophys_ecrit('knf',klev,'Kz final','m2/s',kn(:,1:klev)) +call iophys_ecrit('kqf',klev,'Kz final','m2/s',kq(:,1:klev)) +endif +#endif + + +ENDIF + + +! print*,'OK2' + RETURN + END SUBROUTINE yamada_c