!! !! This module computes diffusion coefficients for continental points. !! !! @author Marie-Alice Foujols and Jan Polcher !! @Version : $Revision: 1.35 $, $Date: 2010/04/07 09:16:40 $ !! !! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/diffuco.f90,v 1.35 2010/04/07 09:16:40 ssipsl Exp $ !! IPSL (2006) !! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC !! MODULE diffuco ! modules used : USE constantes USE constantes_veg USE sechiba_io USE ioipsl USE constantes_co2 USE parallel ! USE WRITE_FIELD_p IMPLICIT NONE ! public routines : ! diffuco_main only PRIVATE PUBLIC :: diffuco_main,diffuco_clear ! ! variables used inside diffuco module : declaration and initialisation ! INTEGER(i_std), PARAMETER :: nlai = 20 LOGICAL, SAVE :: l_first_diffuco = .TRUE. !! Initialisation has to be done one time CHARACTER(LEN=80) :: var_name !! To store variables names for I/O REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: leaf_ci !! intercellular CO2 concentration (ppm) REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rstruct !! architectural resistance REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: raero !! Aerodynamic resistance REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsatt !! Surface saturated humidity !! Nathalie le 28 mars 2006 - sur proposition de Fred Hourdin, ajout !! d'un potentiometre pour regler la resistance de la vegetation REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: rveg_pft ! MM REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: wind !! Wind norm CONTAINS !! !! Main routine for *diffuco* module. !! - called only one time for initialisation !! - called every time step !! - called one more time at last time step for writing _restart_ file !! !! Algorithm: !! - call diffuco_aero for aerodynamic transfer coeficient !! - call diffuco_snow for partial beta coefficient : sublimation !! - call diffuco_inter for partial beta coefficient : interception for each type of vegetation !! - call diffuco_bare for partial beta coefficient : bare soil !! - call diffuco_trans for partial beta coefficient : transpiration for each type of vegetation !! - call diffuco_comb for alpha and beta coefficient !! !! @call diffuco_aero !! @call diffuco_snow !! @call diffuco_inter !! @call diffuco_bare !! @call diffuco_trans !! @call diffuco_comb !! SUBROUTINE diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, u, v, & ! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget ! & zlev, z0, roughheight, temp_sol, temp_air, rau, q_cdrag, qsurf, qair, pb, & & zlev, z0, roughheight, temp_sol, temp_air, rau, q_cdrag, qsurf, qair, q2m, t2m, pb, & & rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, & & swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, & & vbeta , valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id) ! interface description: ! input scalar INTEGER(i_std), INTENT(in) :: kjit !! Time step number INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds LOGICAL, INTENT(in) :: ldrestart_read !! Logical for restart file to read LOGICAL, INTENT(in) :: ldrestart_write !! Logical for restart file to write ! input fields INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0 !! Surface roughness REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective height for roughness REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity ! Ajout Nathalie - declaration q2m REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface level pressure REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rsol !! Bare soil evaporation resistance REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evap_bare_lim !! Beta factor for bare soil evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,... REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,... REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in):: assim_param !! min+max+opt temps, vcmax, vjmax for photosynthesis ! modified fields REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: humrel !! Moisture stress REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag ! Aerodynamic conductance ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: valpha !! Total alpha coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! STOMATE: Beta for CO2 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for transpiration REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Surface resistance for the vegetatuon REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! STOMATE: mean intercellular ci (see enerbil) ! Local ! AJout Nathalie - Juin 2006 !! Beta for fraction of wetted foliage that will transpire REAL(r_std),DIMENSION (kjpindex,nvm) :: vbeta23 INTEGER(i_std) :: ilai CHARACTER(LEN=4) :: laistring ! do initialisation if needed IF (l_first_diffuco) THEN !Config Key = CDRAG_FROM_GCM !Config Desc = Keep cdrag coefficient from gcm. !Config Def = TRUE if q_cdrag on initialization is non zero !Config Help = Set to .TRUE. if you want q_cdrag coming from GCM. !Congig Keep cdrag coefficient from gcm for latent and sensible heat fluxes. IF ( ABS(MAXVAL(q_cdrag)) .LE. EPSILON(q_cdrag)) THEN ldq_cdrag_from_gcm = .FALSE. ELSE ldq_cdrag_from_gcm = .TRUE. ENDIF !MM q_cdrag is always 0 on initialization ?? CALL getin_p('CDRAG_from_GCM', ldq_cdrag_from_gcm) WRITE(numout,*) "ldq_cdrag_from_gcm = ",ldq_cdrag_from_gcm IF (long_print) WRITE (numout,*) ' call diffuco_init ' ! If cdrag is CALL diffuco_init(kjit, ldrestart_read, kjpindex, index, rest_id, q_cdrag) RETURN ENDIF ! ! prepares restart file for the next simulation ! IF (ldrestart_write) THEN IF (long_print) WRITE (numout,*) ' we have to complete restart file with DIFFUCO variables ' var_name= 'rstruct' CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, rstruct, 'scatter', nbp_glo, index_g) var_name= 'raero' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, raero, 'scatter', nbp_glo, index_g) var_name= 'qsatt' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, qsatt, 'scatter', nbp_glo, index_g) ! the following variable is written only if CO2 was calculated IF ( control%ok_co2 ) THEN DO ilai = 1, nlai ! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment... write(laistring,'(i4)') ilai laistring=ADJUSTL(laistring) var_name='leaf_ci_'//laistring(1:LEN_TRIM(laistring)) CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, leaf_ci(:,:,ilai), 'scatter', nbp_glo, index_g) ENDDO ENDIF IF (.NOT.ldq_cdrag_from_gcm) THEN var_name= 'cdrag' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, q_cdrag, 'scatter', nbp_glo, index_g) ENDIF RETURN END IF ! MM wind(:) = SQRT (u(:)*u(:) + v(:)*v(:)) ! ! calculs des differents coefficients ! IF (.NOT.ldq_cdrag_from_gcm) THEN CALL diffuco_aero (kjpindex, kjit, u, v, zlev, z0, roughheight, veget_max, temp_sol, temp_air, & & qsurf, qair, q_cdrag) ENDIF CALL diffuco_raerod (kjpindex, u, v, q_cdrag, raero) ! ! An estimation of the satturated humidity at the surface ! CALL qsatcalc (kjpindex, temp_sol, pb, qsatt) ! ! beta coefficient for sublimation ! CALL diffuco_snow (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, & & snow, frac_nobio, totfrac_nobio, snow_nobio, vbeta1) ! ! beta coefficient for interception ! ! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23 !CALL diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, & ! & qsintveg, qsintmax, rstruct, vbeta2) CALL diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, & & qsintveg, qsintmax, rstruct, vbeta2, vbeta23) ! ! beta coefficient for bare soil ! CALL diffuco_bare (kjpindex, dtradia, u, v, q_cdrag, rsol, evap_bare_lim, evapot, humrel, veget, vbeta4) ! ! beta coefficient for transpiration ! IF ( control%ok_co2 ) THEN ! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget ! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23 !CALL diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, & ! assim_param, ccanopy, & ! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2) CALL diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, q2m, t2m, rau, u, v, q_cdrag, humrel, & assim_param, ccanopy, & veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23) ELSE ! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23 !CALL diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, & ! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2) CALL diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, & veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23) ENDIF ! ! combination of coefficient : alpha and beta coefficient ! ! Ajout qsintmax dans les arguments de la routine.... Nathalie / le 13-03-2006 CALL diffuco_comb (kjpindex, dtradia, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, snow, & & veget, lai, vbeta1, vbeta2, vbeta3, vbeta4, valpha, vbeta, qsintmax) IF ( .NOT. almaoutput ) THEN CALL histwrite(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg) CALL histwrite(hist_id, 'raero', kjit, raero, kjpindex, index) ! Ajouts Nathalie - novembre 2006 CALL histwrite(hist_id, 'cdrag', kjit, q_cdrag, kjpindex, index) CALL histwrite(hist_id, 'Wind', kjit, wind, kjpindex, index) ! Fin ajouts Nathalie !MM CALL histwrite(hist_id, 'qsatt', kjit, qsatt, kjpindex, index) IF ( hist2_id > 0 ) THEN CALL histwrite(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg) CALL histwrite(hist2_id, 'raero', kjit, raero, kjpindex, index) CALL histwrite(hist2_id, 'cdrag', kjit, q_cdrag, kjpindex, index) CALL histwrite(hist2_id, 'Wind', kjit, wind, kjpindex, index) CALL histwrite(hist2_id, 'qsatt', kjit, qsatt, kjpindex, index) ENDIF ELSE ENDIF ! ! IF (long_print) WRITE (numout,*) ' diffuco_main done ' END SUBROUTINE diffuco_main !! Algorithm: !! - dynamic allocation for local array !! SUBROUTINE diffuco_init(kjit, ldrestart_read, kjpindex, index, rest_id, q_cdrag) ! interface description ! input scalar INTEGER(i_std), INTENT (in) :: kjit !! Time step number LOGICAL,INTENT (in) :: ldrestart_read !! Logical for restart file to read INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size INTEGER(i_std),DIMENSION (kjpindex), INTENT (in):: index !! Indeces of the points on the map INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag ! input fields ! output scalar ! output fields ! local declaration INTEGER(i_std) :: ier, jv INTEGER(i_std) :: ilai CHARACTER(LEN=4) :: laistring REAL(r_std),DIMENSION (kjpindex) :: temp ! ! initialisation ! IF (l_first_diffuco) THEN l_first_diffuco=.FALSE. ELSE WRITE (numout,*) ' l_first_diffuco false . we stop ' STOP 'diffuco_init' ENDIF ! allocate only if CO2 is calculated IF ( control%ok_co2 ) THEN ALLOCATE (leaf_ci(kjpindex,nvm,nlai),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in leaf_ci allocation. We stop. We need kjpindex*nvm*nlai words = ',& kjpindex*nvm*nlai STOP 'diffuco_init' END IF ENDIF ALLOCATE (rstruct(kjpindex,nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in rstruct allocation. We stop. We need kjpindex x nvm words = ' ,kjpindex,' x ' ,nvm,& & ' = ',kjpindex*nvm STOP 'diffuco_init' END IF ALLOCATE (raero(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in raero allocation. We stop. We need kjpindex x nvm words = ', kjpindex STOP 'diffuco_init' END IF ALLOCATE (qsatt(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in qsatt allocation. We stop. We need kjpindex x nvm words = ', kjpindex STOP 'diffuco_init' END IF ALLOCATE (wind(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in wind allocation. We stop. We need kjpindex x nvm words = ', kjpindex STOP 'diffuco_init' END IF IF (ldrestart_read) THEN IF (long_print) WRITE (numout,*) ' we have to read a restart file for DIFFUCO variables' var_name='rstruct' CALL ioconf_setatt('UNITS', 's/m') CALL ioconf_setatt('LONG_NAME','Structural resistance') CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., rstruct, "gather", nbp_glo, index_g) ! IF ( MINVAL(rstruct) .EQ. MAXVAL(rstruct) .AND. MAXVAL(rstruct) .EQ. val_exp ) THEN ! DO jv = 1, nvm rstruct(:,jv) = rstruct_const(jv) ENDDO ENDIF var_name='raero' ; CALL ioconf_setatt('UNITS', 's/m') CALL ioconf_setatt('LONG_NAME','Aerodynamic resistance') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN raero(:) = temp(:) ENDIF ENDIF ! var_name='qsatt' ; CALL ioconf_setatt('UNITS', 'g/g') CALL ioconf_setatt('LONG_NAME','Surface saturated humidity') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN qsatt(:) = temp(:) ENDIF ENDIF ! the following variable is read only if CO2 is calculated IF ( control%ok_co2 ) THEN CALL ioconf_setatt('UNITS', 'ppm') CALL ioconf_setatt('LONG_NAME','Leaf CO2') DO ilai = 1, nlai ! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment... write(laistring,'(i4)') ilai laistring=ADJUSTL(laistring) var_name='leaf_ci_'//laistring(1:LEN_TRIM(laistring)) CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE.,leaf_ci(:,:,ilai), "gather", nbp_glo, index_g) ENDDO ! !Config Key = DIFFUCO_LEAFCI !Config Desc = Initial leaf CO2 level if not found in restart !Config Def = 233. !Config Help = The initial value of leaf_ci if its value is not found !Config in the restart file. This should only be used if the model is !Config started without a restart file. ! CALL setvar_p (leaf_ci, val_exp,'DIFFUCO_LEAFCI', 233._r_std) ENDIF ! IF (.NOT.ldq_cdrag_from_gcm) THEN var_name= 'cdrag' CALL ioconf_setatt('LONG_NAME','Drag coefficient for LE and SH') CALL ioconf_setatt('UNITS', '-') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN q_cdrag(:) = temp(:) ENDIF ENDIF ENDIF ENDIF ! ! Ajouts Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin ! ALLOCATE (rveg_pft(nvm),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) & ' error in rveg_pft allocation. We stop. We need nvm words = ', nvm STOP 'diffuco_init' END IF ! !Config Key = RVEG_PFT !Config Desc = Artificial parameter to increase or decrease canopy resistance. !Config Def = 1. !Config Help = This parameter is set by PFT. rveg_pft(:) = 1. CALL getin_p('RVEG_PFT', rveg_pft) WRITE(numout,*) 'DANS DIFFUCO_INIT , RVEG_PFT=',rveg_pft IF (long_print) WRITE (numout,*) ' diffuco_init done ' END SUBROUTINE diffuco_init SUBROUTINE diffuco_clear() l_first_diffuco=.TRUE. IF (ALLOCATED (leaf_ci)) DEALLOCATE (leaf_ci) IF (ALLOCATED (rstruct)) DEALLOCATE (rstruct) IF (ALLOCATED (raero)) DEALLOCATE (raero) IF (ALLOCATED (rveg_pft)) DEALLOCATE (rveg_pft) END SUBROUTINE diffuco_clear !! This routine computes aerothermic coefficient if required !! see logical *ldq_cdrag_from_gcm* !! SUBROUTINE diffuco_aero (kjpindex, kjit, u, v, zlev, z0, roughheight, veget_max, temp_sol, temp_air, & & qsurf, qair, q_cdrag) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex, kjit !! Domain size ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0 !! Surface roughness REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective roughness height REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity ! output scalar ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag ! Aerodynamic conductance ! local declaration INTEGER(i_std) :: ji, jv REAL(r_std) :: speed, zg, zdphi, ztvd, ztvs, zdu2 REAL(r_std) :: zri, cd_neut, zscf, cd_tmp ! initialisation ! test if we have to work with q_cdrag or to calcul it DO ji=1,kjpindex ! ! 1. computes wind speed ! speed = wind(ji) ! ! 2. computes geopotentiel ! zg = zlev(ji) * cte_grav zdphi = zg/cp_air ! ! 3. virtual air temperature at the surface ! ztvd = (temp_air(ji) + zdphi / (un + rvtmp2 * qair(ji))) * (un + retv * qair(ji)) ! ! 4. virtual surface temperature ! ztvs = temp_sol(ji) * (un + retv * qsurf(ji)) ! ! 5. squared wind shear ! zdu2 = MAX(cepdu2,speed**2) ! ! 6. Richardson number ! zri = zg * (ztvd - ztvs) / (zdu2 * ztvd) zri = MAX(MIN(zri,5.),-5.) ! ! 7. Computing the drag coefficient ! We add the add the height of the vegetation to the level height to take into account ! that the level seen by the vegetation is actually the top of the vegetation. Then we ! we can subtract the displacement height. ! cd_neut = (ct_karman / LOG( (zlev(ji) + roughheight(ji)) / z0(ji) )) ** 2 ! ! 7.1 Stable case ! IF (zri .GE. zero) THEN zscf = SQRT(un + cd * ABS(zri)) cd_tmp=cd_neut/(un + trois * cb * zri * zscf) ELSE ! ! 7.2 Unstable case ! zscf = un / (un + trois * cb * cc * cd_neut * SQRT(ABS(zri) * & & ((zlev(ji) + roughheight(ji)) / z0(ji)))) cd_tmp=cd_neut * (un - trois * cb * zri * zscf) ENDIF ! dont let it go to low else the surface uncouples q_cdrag(ji) = MAX(cd_tmp, 1.e-4/MAX(speed,min_wind)) !! !! In some situations it might be useful to give an upper limit on the cdrag as well. !! The line here should then be uncommented. !! q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind)) END DO IF (long_print) WRITE (numout,*) ' not ldqcdrag_from_gcm : diffuco_aero done ' END SUBROUTINE diffuco_aero !! This routine computes partial beta coefficient : snow sublimation !! SUBROUTINE diffuco_snow (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, & & snow, frac_nobio, totfrac_nobio, snow_nobio, vbeta1) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,... REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,... REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio!! Total fraction of ice+lakes+cities+... ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation ! local declaration REAL(r_std) :: subtest, zrapp, speed, vbeta1_add INTEGER(i_std) :: ji, jv ! ! 1. beta coefficient for sublimation for snow on vegetation ! DO ji=1,kjpindex ! Fraction of mesh that can sublimate snow vbeta1(ji) = (un - totfrac_nobio(ji)) * MAX( MIN(snow(ji)/snowcri,un), zero) ! ! -- Limitation of sublimation in case of snow amounts smaller than ! the atmospheric demande. ! speed = MAX(min_wind, wind(ji)) ! subtest = dtradia * vbeta1(ji) * speed * q_cdrag(ji) * rau(ji) * & & ( qsatt(ji) - qair(ji) ) ! IF ( subtest .GT. zero ) THEN zrapp = snow(ji) / subtest IF ( zrapp .LT. un ) THEN vbeta1(ji) = vbeta1(ji) * zrapp ENDIF ENDIF ! END DO ! ! 2. add beta coefficient for other surface types. ! DO jv = 1, nnobio !!$ ! !!$ IF ( jv .EQ. iice ) THEN !!$ ! !!$ ! Land ice is of course a particular case !!$ ! !!$ DO ji=1,kjpindex !!$ vbeta1(ji) = vbeta1(ji) + frac_nobio(ji,jv) !!$ ENDDO !!$ ! !!$ ELSE ! DO ji=1,kjpindex ! vbeta1_add = frac_nobio(ji,jv) * MAX( MIN(snow_nobio(ji,jv)/snowcri,un), zero) ! ! -- Limitation of sublimation in case of snow amounts smaller than ! the atmospheric demand. ! speed = MAX(min_wind, wind(ji)) ! subtest = dtradia * vbeta1_add * speed * q_cdrag(ji) * rau(ji) * & & ( qsatt(ji) - qair(ji) ) ! IF ( subtest .GT. zero ) THEN zrapp = snow_nobio(ji,jv) / subtest IF ( zrapp .LT. un ) THEN vbeta1_add = vbeta1_add * zrapp ENDIF ENDIF ! vbeta1(ji) = vbeta1(ji) + vbeta1_add ! ENDDO !!$ ! !!$ ENDIF ! ENDDO IF (long_print) WRITE (numout,*) ' diffuco_snow done ' END SUBROUTINE diffuco_snow !! This routine computes partial beta coefficient : interception for each type of vegetation !! ! Nathalie - Juin 2006 - Introduction de vbeta23 !SUBROUTINE diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, & ! & qsintveg, qsintmax, rstruct, vbeta2) SUBROUTINE diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, & & qsintveg, qsintmax, rstruct, vbeta2, vbeta23) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! vegetation fraction for each type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! STOMATE: architectural resistance ! output fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss ! AJout Nathalie - Juin 2006 !! Beta for fraction of wetted foliage that will transpire REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23 ! Fin ajout Nathalie ! local declaration INTEGER(i_std) :: ji, jv REAL(r_std) :: zqsvegrap, ziltest, zrapp, speed ! ! Correction Nathalie - Initialisation des vbeta2x vbeta2(:,:) = zero ! Ajout Nathalie - Juin 2006 vbeta23(:,:) = zero ! Fin ajout Nathalie ! DO jv = 1,nvm ! ! 1. beta coefficient for vegetation interception ! DO ji=1,kjpindex IF (veget(ji,jv) .GT. min_sechiba .AND. qsintveg(ji,jv) .GT. zero ) THEN zqsvegrap = zero IF (qsintmax(ji,jv) .GT. min_sechiba ) THEN zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv)) END IF ! speed = MAX(min_wind, wind(ji)) ! -- Interception loss: IL vbeta2(ji,jv) = veget(ji,jv) * zqsvegrap * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv))) ! ! -- Limitation of IL by the water stored on the leaf. ! A first approximation of IL is obtained with the old values of ! qair and qsol_sat: function of temp-sol and pb. (see call of qsatcalc) ! ziltest = dtradia * vbeta2(ji,jv) * speed * q_cdrag(ji) * rau(ji) * & & ( qsatt(ji) - qair(ji) ) IF ( ziltest .GT. zero ) THEN zrapp = qsintveg(ji,jv) / ziltest IF ( zrapp .LT. un ) THEN ! Ajout Nathalie - Juin 2006 vbeta23(ji,jv) = MAX(vbeta2(ji,jv) - vbeta2(ji,jv) * zrapp, 0.) ! Fin ajout Nathalie vbeta2(ji,jv) = vbeta2(ji,jv) * zrapp ENDIF ENDIF END IF END DO END DO IF (long_print) WRITE (numout,*) ' diffuco_inter done ' END SUBROUTINE diffuco_inter !! This routine computes partial beta coefficient : bare soil !! SUBROUTINE diffuco_bare (kjpindex, dtradia, u, v, q_cdrag, rsol, evap_bare_lim, evapot, humrel, veget, vbeta4) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rsol !! resistance for bare soil evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evap_bare_lim !! Beta factor for bare soil evaporation REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Type of vegetation fraction ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation ! local declaration INTEGER(i_std) :: ji, jv REAL(r_std) :: speed IF ( .NOT. control%hydrol_cwrr ) THEN DO ji = 1, kjpindex ! vbeta4(ji) = zero ! ! 1. Soil resistance and beta for bare soil ! note: veget ( ,1) contains the fraction of bare soil ! see hydrol module ! IF (veget(ji,1) .GE. min_sechiba) THEN ! speed = MAX(min_wind, wind(ji)) ! ! Correction Nathalie de Noblet - le 27 Mars 2006 ! Selon recommandation de Frederic Hourdin: supprimer humrel dans formulation vbeta4 !vbeta4(ji) = veget(ji,1) *humrel(ji,1)* (un / (un + speed * q_cdrag(ji) * rsol(ji))) ! Nathalie - le 28 mars 2006 - vbeta4 n'etait pas calcule en fonction de ! rsol mais de rsol_cste * hdry! Dans ce cas inutile de calculer rsol(ji)!! vbeta4(ji) = veget(ji,1) * (un / (un + speed * q_cdrag(ji) * rsol(ji))) ! ENDIF ! END DO ELSE DO ji = 1, kjpindex vbeta4(ji) = evap_bare_lim(ji) END DO ENDIF IF (long_print) WRITE (numout,*) ' diffuco_bare done ' END SUBROUTINE diffuco_bare !! This routine computes partial beta coefficient : transpiration for each type of vegetation !! ! Nathalie - Juin 2006 - introduction de vbeta23 !SUBROUTINE diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, & ! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2) SUBROUTINE diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, & veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Short wave net flux in REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Type of vegetation fraction REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. vegetation fraction (LAI -> infty) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! STOMATE ! AJout Nathalie - Juin 2006 !! Beta for fraction of wetted foliage that will transpire REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 ! Fin ajout Nathalie ! output fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Surface resistance of vegetation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! STOMATE REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! STOMATE ! local declaration INTEGER(i_std) :: ji, jv REAL(r_std) :: speed REAL(r_std), DIMENSION(kjpindex) :: zdefconc, zqsvegrap REAL(r_std), DIMENSION(kjpindex) :: qsatt ! ! 1. Moisture concentration at the leaf level. ! CALL qsatcalc (kjpindex, temp_air, pb, qsatt) zdefconc(:) = rau(:) * MAX( qsatt(:) - qair(:), zero ) ! ! 2. beta coefficient for vegetation transpiration ! DO jv = 1,nvm rveget(:,jv) = undef_sechiba vbeta3(:,jv) = zero zqsvegrap(:) = zero DO ji = 1, kjpindex speed = MAX(min_wind, wind(ji)) IF (qsintmax(ji,jv) .GT. min_sechiba) THEN zqsvegrap(ji) = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv)) ENDIF IF ( ( veget(ji,jv)*lai(ji,jv) .GT. min_sechiba ) .AND. & ( kzero(jv) .GT. min_sechiba ) .AND. & ( swnet(ji) .GT. min_sechiba ) ) THEN rveget(ji,jv) = (( swnet(ji) + rayt_cste ) / swnet(ji) ) & * ((defc_plus + (defc_mult * zdefconc(ji) )) / kzero(jv)) * (un / lai(ji,jv)) ! Corrections Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin ! Introduction d'un potentiometre (rveg_pft) pour regler la somme rveg+rstruct !vbeta3(ji,jv) = veget(ji,jv) * (un - zqsvegrap(ji)) * humrel(ji,jv) * & ! (un / (un + speed * q_cdrag(ji) * (rveget(ji,jv) + rstruct(ji,jv)))) vbeta3(ji,jv) = veget(ji,jv) * (un - zqsvegrap(ji)) * humrel(ji,jv) * & (un / (un + speed * q_cdrag(ji) * (rveg_pft(jv)*(rveget(ji,jv) + rstruct(ji,jv))))) ! Fin ajout Nathalie ! Ajout Nathalie - Juin 2006 vbeta3(ji,jv) = vbeta3(ji,jv) + MIN( vbeta23(ji,jv), & veget(ji,jv) * zqsvegrap(ji) * humrel(ji,jv) * & (un / (un + speed * q_cdrag(ji) * (rveg_pft(jv)*(rveget(ji,jv) + rstruct(ji,jv)))))) ! Fin ajout Nathalie ENDIF ENDDO ENDDO ! STOMATE cimean(:,:) = zero vbetaco2(:,:) = zero IF (long_print) WRITE (numout,*) ' diffuco_trans done ' END SUBROUTINE diffuco_trans !! This routine computes partial beta coefficient : transpiration for each type of vegetation !! STOMATE: this routine now calculates also the assimilation using the Farqhuar & al (1980) formulation !! ! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget ! Nathalie - Juin 2006 - introduction de vbeta23 !SUBROUTINE diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, & ! assim_param, ccanopy, & ! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2) SUBROUTINE diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, q2m, t2m, rau, u, v, q_cdrag, humrel, & assim_param, ccanopy, & veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Downwelling short wave flux REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity ! Ajout Nathalie - Juin 2006 - declaration q2m REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag !! min+max+opt temps, vcmax, vjmax for photosynthesis REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! STOMATE: CO2 concentration inside the canopy REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Type of vegetation fraction REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. vegetation fraction (LAI -> infty) REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation ! AJout Nathalie - Juin 2006 !! Beta for fraction of wetted foliage that will transpire REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 ! Fin ajout Nathalie ! output fields REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Surface resistance of vegetation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! STOMATE REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! STOMATE REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! STOMATE ! local declaration REAL(r_std),DIMENSION (kjpindex,nvm) :: vcmax REAL(r_std),DIMENSION (kjpindex,nvm) :: vjmax REAL(r_std),DIMENSION (kjpindex,nvm) :: tmin REAL(r_std),DIMENSION (kjpindex,nvm) :: topt REAL(r_std),DIMENSION (kjpindex,nvm) :: tmax INTEGER(i_std) :: ji, jv, jl REAL(r_std), DIMENSION(kjpindex) :: leaf_ci_lowest INTEGER(i_std), DIMENSION(kjpindex) :: ilai REAL(r_std), DIMENSION(kjpindex) :: zqsvegrap REAL(r_std) :: speed ! STOMATE: LOGICAL, DIMENSION(kjpindex) :: assimilate, calculate INTEGER(i_std) :: nic,inic,icinic INTEGER(i_std), DIMENSION(kjpindex) :: index_calc INTEGER(i_std) :: nia,inia,nina,inina,iainia INTEGER(i_std), DIMENSION(kjpindex) :: index_assi,index_non_assi REAL(r_std), PARAMETER :: laimax = 12. REAL(r_std), PARAMETER :: xc4_1 = .83 REAL(r_std), PARAMETER :: xc4_2 = .93 REAL(r_std), DIMENSION(kjpindex) :: vc2, vj2 REAL(r_std), DIMENSION(kjpindex) :: assimi REAL(r_std) :: x_1,x_2,x_3,x_4,x_5,x_6 REAL(r_std), DIMENSION(kjpindex) :: gstop, gs REAL(r_std), DIMENSION(kjpindex) :: Kc, Ko, CP REAL(r_std), DIMENSION(kjpindex) :: vc, vj REAL(r_std), DIMENSION(kjpindex) :: kt, rt REAL(r_std), DIMENSION(kjpindex) :: air_relhum REAL(r_std), DIMENSION(kjpindex) :: water_lim, temp_lim REAL(r_std), DIMENSION(kjpindex) :: gstot REAL(r_std), DIMENSION(kjpindex) :: assimtot REAL(r_std), DIMENSION(kjpindex) :: leaf_gs_top !! stomatal conductance at topmost level REAL(r_std), DIMENSION(nlai+1) :: laitab !! tabulated LAI steps REAL(r_std), DIMENSION(kjpindex) :: qsatt REAL(r_std), DIMENSION(nvm,nlai) :: light !! fraction of light that gets through REAL(r_std), DIMENSION(kjpindex) :: ci_gs REAL(r_std) :: cresist !! coefficient for resistances ! ! calculate LAI steps ! DO jl = 1, nlai+1 laitab(jl) = laimax*(EXP(.15*REAL(jl-1,r_std))-1.)/(EXP(.15*REAL(nlai,r_std))-1.) ENDDO ! ! calculate light fraction that comes through at a given LAI for each vegetation type ! DO jl = 1, nlai ! DO jv = 1, nvm ! light(jv,jl) = exp( -ext_coef(jv)*laitab(jl) ) ! ENDDO ! ENDDO ! ! 1. Photosynthesis parameters ! ! ! temperatures in K ! tmin(:,:) = assim_param(:,:,itmin) tmax(:,:) = assim_param(:,:,itmax) topt(:,:) = assim_param(:,:,itopt) ! vcmax(:,:) = assim_param(:,:,ivcmax) vjmax(:,:) = assim_param(:,:,ivjmax) ! ! estimation of relative humidity of the air ! ! correction Nathalie, on utilise q2m/t2m au lieu de qair - Juin 2006 ! CALL qsatcalc (kjpindex, temp_air, pb, qsatt) ! air_relhum(:) = & ! ( qair(:) * pb(:) / (0.622+qair(:)*0.378) ) / & ! ( qsatt(:)*pb(:) / (0.622+qsatt(:)*0.378 ) ) CALL qsatcalc (kjpindex, t2m, pb, qsatt) air_relhum(:) = & ( q2m(:) * pb(:) / (0.622+q2m(:)*0.378) ) / & ( qsatt(:)*pb(:) / (0.622+qsatt(:)*0.378 ) ) ! DO jv = 1,nvm ! ! 2. beta coefficient for vegetation transpiration ! rstruct(:,jv) = rstruct_const(jv) rveget(:,jv) = undef_sechiba ! vbeta3(:,jv) = zero vbetaco2(:,jv) = zero ! cimean(:,jv) = ccanopy(:) ! ! mask that contains points where there is photosynthesis ! nia=0 nina=0 ! DO ji=1,kjpindex ! IF ( ( lai(ji,jv) .GT. 0.01 ) .AND. & ( veget_max(ji,jv) .GT. 1.E-8 ) ) THEN IF ( ( veget(ji,jv) .GT. 1.E-8 ) .AND. & ( swdown(ji) .GT. min_sechiba ) .AND. & ( temp_air(ji) .GT. tmin(ji,jv) ) .AND. & ( temp_air(ji) .LT. tmax(ji,jv) ) .AND. & ( humrel(ji,jv) .GT. min_sechiba ) ) THEN ! assimilate(ji) = .TRUE. nia=nia+1 index_assi(nia)=ji ! ELSE ! assimilate(ji) = .FALSE. nina=nina+1 index_non_assi(nina)=ji ! ENDIF ELSE ! assimilate(ji) = .FALSE. nina=nina+1 index_non_assi(nina)=ji ! ENDIF ! ENDDO ! gstot(:) = zero assimtot(:) = zero ! zqsvegrap(:) = zero WHERE (qsintmax(:,jv) .GT. min_sechiba) zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv)) ENDWHERE ! WHERE ( assimilate(:) ) water_lim(:) = MIN( 2.*humrel(:,jv), 1. ) ENDWHERE ! give a default value of ci for all pixel that do not assimilate DO jl=1,nlai DO inina=1,nina leaf_ci(index_non_assi(inina),jv,jl) = ccanopy(index_non_assi(inina)) * .667_r_std ENDDO ENDDO ! ilai(:) = 1 ! ! Here is calculated photosynthesis (Farqhuar et al. 80) ! and stomatal conductance (Ball & al. 86) ! ! Calculating temperature dependent parameters ! IF ( is_c4(jv) ) THEN ! ! Case of C4 plants ! IF (nia .GT. 0) then !OCL NOVREC DO inia=1,nia ! x_1 = 0.177 * EXP( 0.069*(temp_air(index_assi(inia))-tp_00) ) ! = 2.0**(((temp_air(index_assi(inia))-tp_00)-25.0)/10.0) ! kt(index_assi(inia)) = 0.7 * x_1 * 1.e6 rt(index_assi(inia)) = 0.8 * x_1 / & ( 1.0 + EXP(1.3*(temp_air(index_assi(inia))-tmax(index_assi(inia),jv))) ) ! vc(index_assi(inia)) = vcmax(index_assi(inia),jv) & * 0.39 * x_1 * water_lim(index_assi(inia)) / & ! * 0.39 * x_1 / & ( (1.0+EXP(0.3*(tmin(index_assi(inia),jv)-temp_air(index_assi(inia))))) & * (1.0+EXP(0.3*(temp_air(index_assi(inia))-topt(index_assi(inia),jv)))) ) ! ENDDO ENDIF ! IF (nina .GT. 0) then ! !OCL NOVREC DO inina=1,nina ! kt(index_non_assi(inina)) = zero rt(index_non_assi(inina)) = zero vc(index_non_assi(inina)) = zero ! ENDDO ! ENDIF ! ELSE ! ! Case of C3 plants ! IF (nia .GT. 0) then ! !OCL NOVREC DO inia=1,nia ! temp_lim(index_assi(inia)) = & (temp_air(index_assi(inia))-tmin(index_assi(inia),jv)) * & (temp_air(index_assi(inia))-tmax(index_assi(inia),jv)) temp_lim(index_assi(inia)) = temp_lim(index_assi(inia)) / & (temp_lim(index_assi(inia))-(temp_air(index_assi(inia))-& topt(index_assi(inia),jv))**2) ! Kc(index_assi(inia)) = 39.09 * EXP(.085*(temp_air(index_assi(inia))-tp_00)) ! Ko(index_assi(inia)) = 2.412 * 210000. & * EXP(.085*(temp_air(index_assi(inia))-tmin(index_assi(inia),jv))) / & (temp_air(index_assi(inia))-tmin(index_assi(inia),jv)) ! CP(index_assi(inia)) = 42. * EXP( 9.46*(temp_air(index_assi(inia))-tp_00-25.)/& temp_air(index_assi(inia)) ) ! vc(index_assi(inia)) = vcmax(index_assi(inia),jv) * & temp_lim(index_assi(inia)) * water_lim(index_assi(inia)) ! temp_lim(index_assi(inia)) vj(index_assi(inia)) = vjmax(index_assi(inia),jv) * & temp_lim(index_assi(inia)) * water_lim(index_assi(inia)) ! temp_lim(index_assi(inia)) ! ENDDO ! ENDIF ! IF (nina .GT. 0) then ! !OCL NOVREC DO inina=1,nina ! temp_lim(index_non_assi(inina)) = zero Kc(index_non_assi(inina)) = zero Ko(index_non_assi(inina)) = zero CP(index_non_assi(inina)) = zero ! vc(index_non_assi(inina)) = zero vj(index_non_assi(inina)) = zero ! ENDDO ! ENDIF ! ENDIF ! C3/C4 ! ! estimate assimilation and conductance for each LAI level ! DO jl = 1, nlai ! nic=0 ! calculate(:) = .FALSE. ! IF (nia .GT. 0) then ! !OCL NOVREC DO inia=1,nia ! calculate(index_assi(inia)) = (laitab(jl) .LE. lai(index_assi(inia),jv) ) ! IF ( calculate(index_assi(inia)) ) THEN ! nic=nic+1 index_calc(nic)=index_assi(inia) ! ENDIF ! ENDDO ! ENDIF ! ! Vmax is scaled into the canopy due to reduction of nitrogen ! x_1 = ( un - .7_r_std * ( un - light(jv,jl) ) ) ! IF ( nic .GT. 0 ) THEN ! DO inic=1,nic ! vc2(index_calc(inic)) = vc(index_calc(inic)) * x_1 vj2(index_calc(inic)) = vj(index_calc(inic)) * x_1 ! ENDDO ! ENDIF ! IF ( is_c4(jv) ) THEN ! ! assimilation for C4 plants (Collatz & al. 91) ! DO ji = 1, kjpindex ! assimi(ji) = 0. ! ENDDO ! IF (nic .GT. 0) THEN ! !OCL NOVREC DO inic=1,nic ! x_1 = - ( vc2(index_calc(inic)) + 0.092 * 2.3* swdown(index_calc(inic)) * & ext_coef(jv) * light(jv,jl) ) x_2 = vc2(index_calc(inic)) * 0.092 * 2.3 * swdown(index_calc(inic)) * & ext_coef(jv) * light(jv,jl) x_3 = ( -x_1 - sqrt( x_1*x_1 - 4.0 * xc4_1 * x_2 ) ) / (2.0*xc4_1) x_4 = - ( x_3 + kt(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) * & 1.0e-6 ) x_5 = x_3 * kt(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) * 1.0e-6 assimi(index_calc(inic)) = ( -x_4 - sqrt( x_4*x_4 - 4. * xc4_2 * x_5 ) ) / (2.*xc4_2) assimi(index_calc(inic)) = assimi(index_calc(inic)) - & rt(index_calc(inic)) ! ENDDO ! ENDIF ! ELSE ! ! assimilation for C3 plants (Farqhuar & al. 80) ! DO ji = 1, kjpindex ! assimi(ji) = 0. ! ENDDO ! IF (nic .GT. 0) THEN ! !OCL NOVREC DO inic=1,nic ! x_1 = vc2(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) / & ( leaf_ci(index_calc(inic),jv,jl) + Kc(index_calc(inic)) * & ( 1._r_std + 210000._r_std / Ko(index_calc(inic)) ) ) x_2 = .8855_r_std*swdown(index_calc(inic))*ext_coef(jv)*light(jv,jl) x_3 = x_2+vj2(index_calc(inic)) x_4 = ( x_3 - sqrt( x_3*x_3 - (4._r_std*.7_r_std*x_2*vj2(index_calc(inic))) ) ) / & (2._r_std*.7_r_std) x_5 = x_4 * leaf_ci(index_calc(inic),jv,jl) / & ( 4.5_r_std * leaf_ci(index_calc(inic),jv,jl) + & 10.5_r_std*CP(index_calc(inic)) ) x_6 = MIN( x_1, x_5 ) assimi(index_calc(inic)) = x_6 * ( 1._r_std - CP(index_calc(inic))/& leaf_ci(index_calc(inic),jv,jl) ) - .011_r_std * vc2(index_calc(inic)) ! ENDDO ! ENDIF ! ENDIF ! IF (nic .GT. 0) THEN ! !OCL NOVREC !cdir NODEP DO inic=1,nic ! ! estimate conductance (Ball & al. 86) ! icinic=index_calc(inic) ! gs(icinic) = water_lim(icinic) * & gs(icinic) = & ( gsslope(jv) * assimi(icinic) * & air_relhum(icinic) / ccanopy(icinic) ) & + gsoffset(jv) gs(icinic) = MAX( gs(icinic), gsoffset(jv) ) ENDDO ! DO inic=1,nic icinic=index_calc(inic) ! ! the new ci is calculated with ! dci/dt=(ccanopy-ci)*gs/1.6-A ! ci=ci+((ccanopy(icinic)-ci)*gs/1.6-& ! assimi(icinic))*dtradia ! we verify that ci is not out of possible values ! ci_gs(icinic) = MIN( ccanopy(icinic), MAX( CP(icinic), & ( ccanopy(icinic) - 1.6_r_std * assimi(icinic) / & gs(icinic) ) ) ) - leaf_ci(icinic,jv,jl) ENDDO !cdir NODEP DO inic=1,nic icinic=index_calc(inic) !to avoid some problem of numerical stability, the leaf_ci is bufferized leaf_ci(icinic,jv,jl) = leaf_ci(icinic,jv,jl) + ci_gs(icinic)/6. ENDDO ! DO inic=1,nic icinic=index_calc(inic) ! ! this might be the last level for which Ci is calculated. Store it for ! initialization of the remaining levels of the Ci array. ! leaf_ci_lowest(icinic) = leaf_ci(icinic,jv,jl) ENDDO ! !cdir NODEP DO inic=1,nic icinic=index_calc(inic) ! ! total assimilation and conductance assimtot(icinic) = assimtot(icinic) + & assimi(icinic) * (laitab(jl+1)-laitab(jl)) gstot(icinic) = gstot(icinic) + & gs(icinic) * (laitab(jl+1)-laitab(jl)) ! ilai(icinic) = jl ! ENDDO ! ENDIF ! ! keep stomatal conductance of topmost level ! IF ( jl .EQ. 1 ) THEN ! leaf_gs_top(:) = 0. ! IF ( nic .GT. 0 ) then ! !OCL NOVREC DO inic=1,nic ! leaf_gs_top(index_calc(inic)) = gs(index_calc(inic)) ! ENDDO ! ENDIF ! ENDIF ! IF (nia .GT. 0) THEN ! !OCL NOVREC DO inia=1,nia ! IF ( .NOT. calculate(index_assi(inia)) ) THEN ! ! a) for plants that are doing photosynthesis, but whose LAI is lower ! than the present LAI step, initialize it to the Ci of the lowest ! canopy level ! leaf_ci(index_assi(inia),jv,jl) = leaf_ci_lowest(index_assi(inia)) ! ENDIF ! ENDDO ! ENDIF ! ENDDO ! loop over LAI steps ! ! final calculations: resistances ! IF (nia .GT. 0) THEN ! !OCL NOVREC !cdir NODEP DO inia=1,nia ! iainia=index_assi(inia) ! ! conversion from mmol/m2/s to m/s ! gstot(iainia) = .0244*(temp_air(iainia)/tp_00)*& (1013./pb(iainia))*gstot(iainia) gstop(iainia) = .0244 * (temp_air(iainia)/tp_00)*& (1013./pb(iainia))*leaf_gs_top(iainia)*& laitab(ilai(iainia)+1) ! rveget(iainia,jv) = 1./gstop(iainia) ! ENDDO ! DO inia=1,nia ! iainia=index_assi(inia) ! ! rstruct is the difference between rtot (=1./gstot) and rveget ! ! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif !rstruct(iainia,jv) = 1./gstot(iainia) - & ! rveget(iainia,jv) rstruct(iainia,jv) = MAX( 1./gstot(iainia) - & rveget(iainia,jv), min_sechiba) ! ENDDO ! DO inia=1,nia ! iainia=index_assi(inia) ! speed = MAX(min_wind, wind(index_assi(inia))) ! ! beta for transpiration ! ! Corrections Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin ! Introduction d'un potentiometre (rveg_pft) pour regler la somme rveg+rstruct !vbeta3(iainia,jv) = veget_max(iainia,jv) * & ! (un - zqsvegrap(iainia)) * & ! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + & ! rstruct(iainia,jv)))) cresist=(un / (un + speed * q_cdrag(iainia) * & (rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv))))) vbeta3(iainia,jv) = veget_max(iainia,jv) * & (un - zqsvegrap(iainia)) * cresist + & !!$ ! Ajout Nathalie - Juin 2006 !!$ vbeta3(iainia,jv) = vbeta3(iainia,jv) + & ! Corrections Nathalie - le 09 novembre 2009 : veget => veget_max ! MIN( vbeta23(iainia,jv), veget(iainia,jv) * & MIN( vbeta23(iainia,jv), veget_max(iainia,jv) * & ! zqsvegrap(iainia) * humrel(iainia,jv) * & zqsvegrap(iainia) * cresist ) ! Fin ajout Nathalie ! ! beta for assimilation. The difference is that surface ! covered by rain (un - zqsvegrap(iainia)) is not taken into account ! 1.6 is conversion for H2O to CO2 conductance ! vbetaco2(iainia,jv) = veget_max(iainia,jv) * & ! (un / (un + q_cdrag(iainia) * & ! (rveget(iainia,jv))))/1.6 ! vbetaco2(iainia,jv) = veget_max(iainia,jv) * & (un / (un + speed * q_cdrag(iainia) * & (rveget(iainia,jv) + rstruct(iainia,jv)))) / 1.6 ! ! cimean is the "mean ci" calculated in such a way that assimilation ! calculated in enerbil is equivalent to assimtot ! cimean(iainia,jv) = ccanopy(iainia) - & assimtot(iainia) / & ( vbetaco2(iainia,jv)/veget_max(iainia,jv) * & rau(iainia) * speed * q_cdrag(iainia)) ! ENDDO ! ENDIF ! END DO ! loop over vegetation types ! IF (long_print) WRITE (numout,*) ' diffuco_trans_co2 done ' END SUBROUTINE diffuco_trans_co2 !! This routine combines previous partial beta coeeficient and calculates !! alpha and complete beta coefficient !! ! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006 SUBROUTINE diffuco_comb (kjpindex, dtradia, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, & & snow, veget, lai, vbeta1, vbeta2, vbeta3 , vbeta4, valpha, vbeta, qsintmax) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! lower air temperature in Kelvin REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index ! Ajout Nathalie / le 13-03-2006 REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation ! modified fields REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta1 !! Beta for sublimation REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta4 !! Beta for Bare soil evaporation REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Soil moisture stress REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2 !! Beta for Interception for REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3 !! Beta for Transpiration ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: valpha !! TotalAlpha coefficient REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient ! local declaration INTEGER(i_std) :: ji, jv REAL(r_std) :: zevtest, zsoil_moist, zrapp REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum REAL(r_std), DIMENSION(kjpindex) :: vegetsum, vegetsum2 REAL(r_std), DIMENSION(kjpindex) :: qsatt LOGICAL, DIMENSION(kjpindex) :: toveg, tosnow REAL(r_std) :: coeff_dew_veg vbeta2sum(:) = 0. vbeta3sum(:) = 0. DO jv = 1, nvm vbeta2sum(:) = vbeta2sum(:) + vbeta2(:,jv) vbeta3sum(:) = vbeta3sum(:) + vbeta3(:,jv) ENDDO ! ! 1. The beta and alpha coefficients are calculated. ! vbeta(:) = un valpha(:) = un ! ! 2. if snow is lower than critical value ! DO ji = 1, kjpindex IF (snow(ji) .LT. snowcri) THEN vbeta(ji) = vbeta4(ji) + vbeta2sum(ji) + vbeta3sum(ji) IF (vbeta(ji) .LT. min_sechiba) THEN vbeta(ji) = zero END IF END IF ENDDO ! ! 3. If we are in presence of dew. ! ! for vectorization: some arrays vegetsum(:) = 0. DO jv = 1, nvm vegetsum(:) = vegetsum(:) + veget(:,jv) ENDDO vegetsum2(:) = 0. DO jv = 2, nvm vegetsum2(:) = vegetsum2(:) + veget(:,jv) ENDDO CALL qsatcalc (kjpindex, temp_sol, pb, qsatt) ! ! 3.1 decide where the water goes (soil, vegetation, or snow) ! when air moisture exceeds saturation. ! toveg(:) = .FALSE. tosnow(:) = .FALSE. DO ji = 1, kjpindex IF ( qsatt(ji) .LT. qair(ji) ) THEN IF (temp_air(ji) .GT. tp_00) THEN ! ! 3.1.1 If it is not freezing dew is put into the ! interception reservoir and on the bare soil. toveg(ji) = .TRUE. ELSE ! ! 3.1.2 If it is freezing water is put into the ! snow reservoir. tosnow(ji) = .TRUE. ENDIF ENDIF END DO ! 3.1.3 now modify valpha and vbetas where necessary. ! ! 3.1.3.1 Soil and snow (2d) ! DO ji = 1, kjpindex IF ( toveg(ji) ) THEN vbeta1(ji) = zero vbeta4(ji) = veget(ji,1) ! Correction Nathalie - le 13-03-2006: le vbeta ne sera calcule qu'une fois tous les vbeta2 redefinis !vbeta(ji) = vegetsum(ji) vbeta(ji) = vbeta4(ji) valpha(ji) = un ENDIF IF ( tosnow(ji) ) THEN vbeta1(ji) = un vbeta4(ji) = zero vbeta(ji) = un valpha(ji) = un ENDIF ENDDO ! ! 3.1.3.2 vegetation (3d) ! DO jv = 1, nvm ! DO ji = 1, kjpindex ! ! Correction Nathalie - 13-03-2006 / si qsintmax=0, vbeta2=0 IF ( toveg(ji) ) THEN IF (qsintmax(ji,jv) .GT. min_sechiba) THEN !MM ! Compute part of dew that can be intercepted by leafs. IF ( lai(ji,jv) .GT. min_sechiba) THEN IF (lai(ji,jv) .GT. 1.5) THEN coeff_dew_veg= & & 0.000017*lai(ji,jv)**5 & & - 0.000824*lai(ji,jv)**4 & & + 0.014843*lai(ji,jv)**3 & & - 0.110112*lai(ji,jv)**2 & & + 0.205673*lai(ji,jv) & & + 0.887773 ELSE coeff_dew_veg=1. ENDIF ELSE coeff_dew_veg=zero ENDIF vbeta2(ji,jv) = coeff_dew_veg*veget(ji,jv) ! vbeta2(ji,jv) = veget(ji,jv) ELSE vbeta2(ji,jv) = zero ENDIF vbeta(ji) = vbeta(ji) + vbeta2(ji,jv) ENDIF IF ( tosnow(ji) ) vbeta2(ji,jv) = zero ! ENDDO ! ENDDO ! ! 3.2 In any case there is no transpiration when air moisture is too high. ! DO jv = 1, nvm DO ji = 1, kjpindex IF ( qsatt(ji) .LT. qair(ji) ) THEN vbeta3(ji,jv) = zero humrel(ji,jv) = zero ENDIF ENDDO ENDDO ! ! 3.2_bis In any case there is no interception loss on bare soil. ! DO ji = 1, kjpindex IF ( qsatt(ji) .LT. qair(ji) ) THEN vbeta2(ji,1) = zero ENDIF ENDDO IF (long_print) WRITE (numout,*) ' diffuco_comb done ' END SUBROUTINE diffuco_comb SUBROUTINE diffuco_raerod (kjpindex, u, v, q_cdrag, raero) ! ! Simply computes the aerodynamic resistance. For the moment it is ! only used as a diagnostic but that may change ! ! IMPLICIT NONE ! INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag ! output filed REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: raero !! Aerodynamic resistance ! ! local declaration INTEGER(i_std) :: ji REAL(r_std) :: speed ! DO ji=1,kjpindex ! speed = MAX(min_wind, wind(ji)) raero(ji) = un / (q_cdrag(ji)*speed) ! ENDDO ! END SUBROUTINE diffuco_raerod END MODULE diffuco