!! !! This module computes soil thermodynamic !! !! @author Marie-Alice Foujols and Jan Polcher !! @Version : $Revision: 45 $, $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $ !! !< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/thermosoil.f90 $ !< $Date: 2011-01-01 21:30:44 +0100 (Sat, 01 Jan 2011) $ !< $Author: mmaipsl $ !< $Revision: 45 $ !! IPSL (2006) !! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC !! MODULE thermosoil ! routines called : restput, restget ! USE ioipsl ! ! modules used : USE constantes USE sechiba_io USE grid USE parallel IMPLICIT NONE ! public routines : ! thermosoil_main PRIVATE PUBLIC :: thermosoil_main,thermosoil_clear ! ! variables used inside thermosoil module : declaration and initialisation ! LOGICAL, SAVE :: l_first_thermosoil=.TRUE. !! Initialisation has to be done one time CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for I/O REAL(r_std), SAVE :: lambda, cstgrnd, lskin, fz1, zalph ! two dimensions array allocated, computed, saved and got in thermosoil module REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! Different levels soil temperature ! one dimension array allocated, computed and used in thermosoil module exclusively REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz !! REAL(r_std), SAVE, DIMENSION (ngrnd) :: zz_coef REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz1 !! REAL(r_std), SAVE, DIMENSION (ngrnd) :: dz2 !! ! one dimension array allocated, computed and used in thermosoil module exclusively REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z1 !! ! two dimensions arrays allocated, computed and used in thermosoil module exclusively REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz1 !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: zdz2 !! ! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! Capacity used for energy_incr calculation REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! Temperature at the beginning of the time step REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at the beginning of the timestep REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow cold content !! REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: wetdiag !! Soil weetness on the thermodynamical levels !! CONTAINS !! Main routine for *thermosoil* module. !! !! Algorithm: !! - call thermosoil_var_init to initialise variables !! - call thermosoil_coef for coefficient !! - call thermosoil_profile for soil profiling !! !! @call thermosoil_var_init !! @call thermosoil_coef !! @call thermosoil_profile !! SUBROUTINE thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexgrnd, & & temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, 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,hist_id !! _Restart_ file and 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*ngrnd), INTENT (in) :: indexgrnd !! Indeces of the points on the 3D map REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow quantity REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostic of relative humidity ! output fields REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! Soil capacity REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout):: stempdiag !! diagnostic temp profile REAL(r_std),DIMENSION (kjpindex,ngrnd) :: temp REAL(r_std),DIMENSION (kjpindex,ngrnd-1) :: temp1 REAL(r_std),DIMENSION (kjpindex) :: temp2 ! ! do initialisation ! IF (l_first_thermosoil) THEN IF (long_print) WRITE (numout,*) ' l_first_thermosoil : call thermosoil_init ' ! ! do needed allocation ! CALL thermosoil_init (kjit, ldrestart_read, kjpindex, index, rest_id) ! ! computes some physical constantes and array depending soil level depth ! CALL thermosoil_var_init (kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag) ! ! computes cgrd and dgrd coefficient from previous time step (restart) ! CALL thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,& & zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa) CALL thermosoil_energy (kjpindex, temp_sol_new, soilcap, .TRUE.) IF (ldrestart_read) THEN IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables' var_name= 'cgrnd' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Cgrnd coefficient.') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g) IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN cgrnd(:,:)=temp1(:,:) ENDIF ENDIF var_name= 'dgrnd' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Dgrnd coefficient.') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g) IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN dgrnd(:,:)=temp1(:,:) ENDIF ENDIF var_name= 'z1' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','?.') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g) IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN z1(:)=temp2(:) ENDIF ENDIF var_name= 'pcapa' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','?.') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN pcapa(:,:)=temp(:,:) ENDIF ENDIF var_name= 'pcapa_en' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','?.') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN pcapa_en(:,:)=temp(:,:) ENDIF ENDIF var_name= 'pkappa' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','?.') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN pkappa(:,:)=temp(:,:) ENDIF ENDIF var_name= 'zdz1' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','?.') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, .TRUE., temp1, "gather", nbp_glo, index_g) IF (MINVAL(temp1) < MAXVAL(temp1) .OR. MAXVAL(temp1) .NE. val_exp) THEN zdz1(:,:)=temp1(:,:) ENDIF ENDIF var_name= 'zdz2' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','?.') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g) IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN zdz2(:,:)=temp(:,:) ENDIF ENDIF var_name='temp_sol_beg' CALL ioconf_setatt('UNITS', 'K') CALL ioconf_setatt('LONG_NAME','Old Surface temperature') IF ( ok_var(var_name) ) THEN CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp2, "gather", nbp_glo, index_g) IF (MINVAL(temp2) < MAXVAL(temp2) .OR. MAXVAL(temp2) .NE. val_exp) THEN temp_sol_beg(:) = temp2(:) ENDIF ENDIF ENDIF 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 THERMOSOIL variables' var_name= 'ptn' CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g) var_name= 'cgrnd' CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g) var_name= 'dgrnd' CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g) var_name= 'z1' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, z1, 'scatter', nbp_glo, index_g) var_name= 'pcapa' CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pcapa, 'scatter', nbp_glo, index_g) var_name= 'pcapa_en' CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pcapa_en, 'scatter', nbp_glo, index_g) var_name= 'pkappa' CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, pkappa, 'scatter', nbp_glo, index_g) var_name= 'zdz1' CALL restput_p(rest_id, var_name, nbp_glo, ngrnd-1, 1, kjit, zdz1, 'scatter', nbp_glo, index_g) var_name= 'zdz2' CALL restput_p(rest_id, var_name, nbp_glo, ngrnd, 1, kjit, zdz2, 'scatter', nbp_glo, index_g) var_name= 'temp_sol_beg' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_beg, 'scatter', nbp_glo, index_g) var_name= 'soilcap' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g) var_name= 'soilflx' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g) ! read in enerbil var_name= 'temp_sol_new' CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, temp_sol_new, 'scatter', nbp_glo, index_g) RETURN END IF ! ! Put the soil wetnesss diagnostic on the levels of the soil temprature ! CALL thermosoil_humlev(kjpindex, shumdiag) ! ! computes profile with previous cgrd and dgrd coefficient ! CALL thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag) CALL thermosoil_energy (kjpindex, temp_sol_new, soilcap, .FALSE.) ! IF ( .NOT. almaoutput ) THEN CALL histwrite(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd) ELSE CALL histwrite(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd) CALL histwrite(hist_id, 'Qg', kjit, soilflx, kjpindex, index) CALL histwrite(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index) CALL histwrite(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index) ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN CALL histwrite(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd) ELSE CALL histwrite(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd) CALL histwrite(hist2_id, 'Qg', kjit, soilflx, kjpindex, index) CALL histwrite(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index) CALL histwrite(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index) ENDIF ENDIF ! ! computes cgrd and dgrd coefficient ! CALL thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,& & zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa) IF (long_print) WRITE (numout,*) ' thermosoil_main done ' END SUBROUTINE thermosoil_main !! Initialisation for *thermosoil* module. !! - does dynamic allocation for local arrays !! - reads _restart_ file or set initial values to a raisonable value !! SUBROUTINE thermosoil_init(kjit, ldrestart_read, kjpindex, index, rest_id) ! 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 ! input fields ! output fields ! local declaration INTEGER(i_std) :: ier ! initialisation IF (l_first_thermosoil) THEN l_first_thermosoil=.FALSE. ELSE WRITE (numout,*) ' l_first_thermosoil false . we stop ' STOP 'thermosoil_init' ENDIF ! two dimensions array allocation ALLOCATE (ptn(kjpindex,ngrnd),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in ptn allocation. We stop. We need ',kjpindex,' fois ',ngrnd,' words = '& & , kjpindex*ngrnd STOP 'thermosoil_init' END IF ! one dimension array allocation ALLOCATE (z1(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in z1 allocation. We STOP. We need ',kjpindex,' words ' STOP 'thermosoil_init' END IF ! two dimension array allocation ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in cgrnd allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '& & , kjpindex*(ngrnd-1) STOP 'thermosoil_init' END IF ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in dgrnd allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '& & , kjpindex*(ngrnd-1) STOP 'thermosoil_init' END IF ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in pcapa allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& & , kjpindex*ngrnd STOP 'thermosoil_init' END IF ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in pkappa allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& & , kjpindex*ngrnd STOP 'thermosoil_init' END IF ALLOCATE (zdz1(kjpindex,ngrnd-1),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in zdz1 allocation. We STOP. We need ',kjpindex,' fois ',ngrnd-1 ,' words = '& & , kjpindex*(ngrnd-1) STOP 'thermosoil_init' END IF ALLOCATE (zdz2(kjpindex,ngrnd),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in zdz2 allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& & , kjpindex*ngrnd STOP 'thermosoil_init' END IF ALLOCATE (surfheat_incr(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in surfheat_incr allocation. We STOP. We need ',kjpindex,' words = '& & , kjpindex STOP 'thermosoil_init' END IF ALLOCATE (coldcont_incr(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in coldcont_incr allocation. We STOP. We need ',kjpindex,' words = '& & , kjpindex STOP 'thermosoil_init' END IF ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in pcapa_en allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& & , kjpindex*ngrnd STOP 'thermosoil_init' END IF ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in ptn_beg allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& & , kjpindex*ngrnd STOP 'thermosoil_init' END IF ALLOCATE (temp_sol_beg(kjpindex),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in temp_sol_beg allocation. We STOP. We need ',kjpindex,' words = '& & , kjpindex STOP 'thermosoil_init' END IF ALLOCATE (wetdiag(kjpindex,ngrnd),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in wetdiag allocation. We STOP. We need ',kjpindex,' fois ',ngrnd ,' words = '& & , kjpindex*ngrnd STOP 'thermosoil_init' END IF ! ! open restart input file done by sechiba_init ! and read data from restart input file for THERMOSOIL process IF (ldrestart_read) THEN IF (long_print) WRITE (numout,*) ' we have to READ a restart file for THERMOSOIL variables' var_name= 'ptn' CALL ioconf_setatt('UNITS', 'K') CALL ioconf_setatt('LONG_NAME','Soil Temperature profile') CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g) ! ! change restart If they were not found in the restart file ! !Config Key = THERMOSOIL_TPRO !Config Desc = Initial soil temperature profile if not found in restart !Config Def = 280. !Config Help = The initial value of the temperature profile in the soil if !Config its value is not found in the restart file. This should only !Config be used if the model is started without a restart file. Here !Config we only require one value as we will assume a constant !Config throughout the column. ! CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std) ENDIF IF (long_print) WRITE (numout,*) ' thermosoil_init done ' END SUBROUTINE thermosoil_init !! Function for distributing the levels !! SUBROUTINE thermosoil_clear() l_first_thermosoil=.TRUE. IF ( ALLOCATED (ptn)) DEALLOCATE (ptn) IF ( ALLOCATED (z1)) DEALLOCATE (z1) IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd) IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd) IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa) IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa) IF ( ALLOCATED (zdz1)) DEALLOCATE (zdz1) IF ( ALLOCATED (zdz2)) DEALLOCATE (zdz2) IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en) IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg) IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg) IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr) IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr) IF ( ALLOCATED (wetdiag)) DEALLOCATE (wetdiag) END SUBROUTINE thermosoil_clear !! !! FUNCTION fz(rk) RESULT (fz_result) ! interface ! input value REAL(r_std), INTENT(in) :: rk ! output value REAL(r_std) :: fz_result fz_result = fz1 * (zalph ** rk - un) / (zalph - un) END FUNCTION fz !! Thermosoil's variables initialisation !! SUBROUTINE thermosoil_var_init(kjpindex, zz, zz_coef, dz1, dz2, pkappa, pcapa, pcapa_en, shumdiag, stempdiag) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostic humidity profile ! output fields REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz !! REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: zz_coef REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz1 !! REAL(r_std), DIMENSION (ngrnd), INTENT(out) :: dz2 !! tailles des couches REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa !! REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pcapa_en REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: pkappa !! REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostic temp profile ! local declaration INTEGER(i_std) :: ier, ji, jg REAL(r_std) :: sum ! ! 0. initialisation ! cstgrnd=SQRT(one_day / pi) lskin = SQRT(so_cond / so_capa * one_day / pi) fz1 = 0.3_r_std * cstgrnd zalph = deux ! ! 1. Computing the depth of the Temperature level, using a ! non dimentional variable x = z/lskin, lskin beeing ! the skin depth ! DO jg=1,ngrnd !!! This needs to be solved soon. Either we allow CPP options in SECHIBA or the VPP !!! fixes its compiler ! !!!#ifdef VPP5000 dz2(jg) = fz(REAL(jg,r_std)-undemi+undemi) - fz(REAL(jg-1,r_std)-undemi+undemi) !!!#else !!! dz2(jg) = fz(REAL(jg,r_std)) - fz(REAL(jg-1,r_std)) !!!#endif ENDDO ! ! 1.2 The undimentional depth is computed. ! ------------------------------------ DO jg=1,ngrnd zz(jg) = fz(REAL(jg,r_std) - undemi) zz_coef(jg) = fz(REAL(jg,r_std)-undemi+undemi) ENDDO ! ! 1.3 Converting to meters. ! -------------------- DO jg=1,ngrnd zz(jg) = zz(jg) / cstgrnd * lskin zz_coef(jg) = zz_coef(jg) / cstgrnd * lskin dz2(jg) = dz2(jg) / cstgrnd * lskin ENDDO ! ! 1.4 Computing some usefull constants. ! -------------------------------- DO jg=1,ngrnd-1 dz1(jg) = un / (zz(jg+1) - zz(jg)) ENDDO lambda = zz(1) * dz1(1) ! ! 1.5 Get the wetness profice on this grid ! ------------------------------------ ! CALL thermosoil_humlev(kjpindex, shumdiag) ! ! 1.6 Thermal conductivity at all levels. ! ---------------------------------- DO jg = 1,ngrnd DO ji = 1,kjpindex pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry) pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry) pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry) ENDDO ENDDO ! ! 2. Diagnostics. ! ----------- WRITE (numout,*) 'diagnostic des niveaux dans le sol' WRITE (numout,*) 'niveaux intermediaires et pleins' sum = zero DO jg=1,ngrnd sum = sum + dz2(jg) WRITE (numout,*) zz(jg),sum ENDDO ! ! 3. Compute the diagnostic temperature profile ! CALL thermosoil_diaglev(kjpindex, stempdiag) ! ! IF (long_print) WRITE (numout,*) ' thermosoil_var_init done ' END SUBROUTINE thermosoil_var_init !! Computation of cgrnd and dgrnd coefficient at the t time-step. !! !! Needs previous time step values. !! SUBROUTINE thermosoil_coef (kjpindex, dtradia, temp_sol_new, snow, ptn, soilcap, soilflx, zz, dz1, dz2, z1, zdz1,& & zdz2, cgrnd, dgrnd, pcapa, pcapa_en, pkappa) ! 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) :: temp_sol_new !! REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: zz !! REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz1 !! REAL(r_std), DIMENSION (ngrnd), INTENT(in) :: dz2 !! REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (in) :: ptn ! output fields REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: z1 !! REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa !! REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pcapa_en !! REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(inout) :: pkappa !! REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !! REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !! REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: zdz1 !! REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(out) :: zdz2 !! ! local declaration INTEGER(i_std) :: ji, jg REAL(r_std), DIMENSION(kjpindex) :: snow_h, zx1, zx2 ! ! objet: computation of cgrnd and dgrnd coefficient at the t time-step. ! ------ ! ! --------------------------------------------------------------- ! Computation of the Cgrd and Dgrd coefficient for the next step: ! --------------------------------------------------------------- ! DO ji = 1,kjpindex snow_h(ji) = snow(ji) / sn_dens ! ! Traitement special pour le premiere couche ! IF ( snow_h(ji) .GT. zz_coef(1) ) THEN pcapa(ji,1) = sn_capa pcapa_en(ji,1) = sn_capa pkappa(ji,1) = sn_cond ELSE IF ( snow_h(ji) .GT. zero ) THEN pcapa_en(ji,1) = sn_capa zx1(ji) = snow_h(ji) / zz_coef(1) zx2(ji) = ( zz_coef(1) - snow_h(ji)) / zz_coef(1) pcapa(ji,1) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet pkappa(ji,1) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet ) ELSE pcapa(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry) pkappa(ji,1) = so_cond_dry + wetdiag(ji,1)*(so_cond_wet - so_cond_dry) pcapa_en(ji,1) = so_capa_dry + wetdiag(ji,1)*(so_capa_wet - so_capa_dry) ENDIF ! DO jg = 2, ngrnd - 2 IF ( snow_h(ji) .GT. zz_coef(jg) ) THEN pcapa(ji,jg) = sn_capa pkappa(ji,jg) = sn_cond pcapa_en(ji,jg) = sn_capa ELSE IF ( snow_h(ji) .GT. zz_coef(jg-1) ) THEN zx1(ji) = (snow_h(ji) - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1)) zx2(ji) = ( zz_coef(jg) - snow_h(ji)) / (zz_coef(jg) - zz_coef(jg-1)) pcapa(ji,jg) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet pkappa(ji,jg) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet ) pcapa_en(ji,jg) = sn_capa ELSE pcapa(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry) pkappa(ji,jg) = so_cond_dry + wetdiag(ji,jg)*(so_cond_wet - so_cond_dry) pcapa_en(ji,jg) = so_capa_dry + wetdiag(ji,jg)*(so_capa_wet - so_capa_dry) ENDIF ENDDO ! ! ENDDO ! DO jg=1,ngrnd DO ji=1,kjpindex zdz2(ji,jg)=pcapa(ji,jg) * dz2(jg)/dtradia ENDDO ENDDO ! DO jg=1,ngrnd-1 DO ji=1,kjpindex zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg) ENDDO ENDDO ! DO ji = 1,kjpindex z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1) cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji) dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji) ENDDO DO jg = ngrnd-1,2,-1 DO ji = 1,kjpindex z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg))) cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji) dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji) ENDDO ENDDO ! --------------------------------------------------------- ! computation of the surface diffusive flux from ground and ! calorific capacity of the ground: ! --------------------------------------------------------- DO ji = 1,kjpindex soilflx(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1)) soilcap(ji) = (zdz2(ji,1) * dtradia + dtradia * (un - dgrnd(ji,1)) * zdz1(ji,1)) z1(ji) = lambda * (un - dgrnd(ji,1)) + un soilcap(ji) = soilcap(ji) / z1(ji) soilflx(ji) = soilflx(ji) + & & soilcap(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dtradia ENDDO IF (long_print) WRITE (numout,*) ' thermosoil_coef done ' END SUBROUTINE thermosoil_coef !! Computation of : the ground temperature evolution !! SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, ptn, stempdiag) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature ! modified fields REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (inout) :: ptn !! Different levels soil temperature ! output fields REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc profile ! local declaration INTEGER(i_std) :: ji, jg ! ! objet: computation of : the ground temperature evolution ! ------ ! ! Method: implicit time integration ! ------- ! Consecutives ground temperatures are related by: ! T(k+1) = C(k) + D(k)*T(k) (1) ! the coefficients C and D are computed at the t-dt time-step. ! Routine structure: ! -new temperatures are computed using (1) ! ! ! surface temperature DO ji = 1,kjpindex ptn(ji,1) = (lambda * cgrnd(ji,1) + temp_sol_new(ji)) / (lambda * (un - dgrnd(ji,1)) + un) ENDDO ! other temperatures DO jg = 1,ngrnd-1 DO ji = 1,kjpindex ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg) ENDDO ENDDO CALL thermosoil_diaglev(kjpindex, stempdiag) IF (long_print) WRITE (numout,*) ' thermosoil_profile done ' END SUBROUTINE thermosoil_profile !! !! !! Diagnostic soil temperature profile on new levels !! !! SUBROUTINE thermosoil_diaglev(kjpindex, stempdiag) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields ! ! modified fields ! ! output fields REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: stempdiag !! Diagnostoc profile ! ! local variable ! INTEGER(i_std) :: ji, jd, jg REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfact ! LOGICAL, PARAMETER :: check=.FALSE. ! ! IF ( .NOT. ALLOCATED(intfact)) THEN ! ALLOCATE(intfact(nbdl, ngrnd)) ! prev_diag = zero DO jd = 1, nbdl lev_diag = diaglev(jd) prev_prog = zero DO jg = 1, ngrnd IF ( jg == ngrnd .AND. (prev_prog + dz2(jg)) < lev_diag ) THEN !! Just make sure we cover the deepest layers lev_prog = lev_diag ELSE lev_prog = prev_prog + dz2(jg) ENDIF intfact(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag) prev_prog = lev_prog ENDDO prev_diag = lev_diag ENDDO ! IF ( check ) THEN WRITE(numout,*) 'thermosoil_diagev -- thermosoil_diaglev -- thermosoil_diaglev --' DO jd = 1, nbdl WRITE(numout,*) jd, '-', intfact(jd,1:ngrnd) ENDDO WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM" DO jd = 1, nbdl WRITE(numout,*) jd, '-', SUM(intfact(jd,1:ngrnd)) ENDDO WRITE(numout,*) 'thermosoil_diaglev -- thermosoil_diaglev -- thermosoil_diaglev --' ENDIF ! ENDIF stempdiag(:,:) = zero DO jg = 1, ngrnd DO jd = 1, nbdl DO ji = 1, kjpindex stempdiag(ji,jd) = stempdiag(ji,jd) + ptn(ji,jg)*intfact(jd,jg) ENDDO ENDDO ENDDO END SUBROUTINE thermosoil_diaglev !! !! !! Put soil wetness on the temperature levels !! !! SUBROUTINE thermosoil_humlev(kjpindex, shumdiag) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size ! input fields ! ! modified fields ! ! output fields REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Diagnostoc profile ! ! local variable ! INTEGER(i_std) :: ji, jd, jg REAL(r_std) :: lev_diag, prev_diag, lev_prog, prev_prog REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: intfactw ! LOGICAL, PARAMETER :: check=.FALSE. ! ! IF ( .NOT. ALLOCATED(intfactw)) THEN ! ALLOCATE(intfactw(ngrnd, nbdl)) ! prev_diag = zero DO jd = 1, ngrnd lev_diag = prev_diag + dz2(jd) prev_prog = zero DO jg = 1, nbdl IF ( jg == nbdl .AND. diaglev(jg) < lev_diag ) THEN !! Just make sure we cover the deepest layers lev_prog = lev_diag ELSE lev_prog = diaglev(jg) ENDIF intfactw(jd,jg) = MAX(MIN(lev_diag,lev_prog)-MAX(prev_diag, prev_prog), zero)/(lev_diag-prev_diag) prev_prog = lev_prog ENDDO prev_diag = lev_diag ENDDO ! IF ( check ) THEN WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --' DO jd = 1, ngrnd WRITE(numout,*) jd, '-', intfactw(jd,1:nbdl) ENDDO WRITE(numout,*) "SUM -- SUM -- SUM SUM -- SUM -- SUM" DO jd = 1, ngrnd WRITE(numout,*) jd, '-', SUM(intfactw(jd,1:nbdl)) ENDDO WRITE(numout,*) 'thermosoil_humlev -- thermosoil_humlev -- thermosoil_humlev --' ENDIF ! ENDIF wetdiag(:,:) = zero DO jg = 1, nbdl DO jd = 1, ngrnd DO ji = 1, kjpindex wetdiag(ji,jd) = wetdiag(ji,jd) + shumdiag(ji,jg)*intfactw(jd,jg) ENDDO ENDDO ENDDO END SUBROUTINE thermosoil_humlev SUBROUTINE thermosoil_energy(kjpindex, temp_sol_new, soilcap, first_call) ! interface description ! input scalar INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size LOGICAL, INTENT (in) :: first_call !! ! input fields ! REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity ! ! modified fields ! ! output fields ! ! local variable ! INTEGER(i_std) :: ji, jg ! ! IF (first_call) THEN DO ji = 1, kjpindex surfheat_incr(ji) = zero coldcont_incr(ji) = zero temp_sol_beg(ji) = temp_sol_new(ji) ! DO jg = 1, ngrnd ptn_beg(ji,jg) = ptn(ji,jg) ENDDO ! ENDDO RETURN ENDIF DO ji = 1, kjpindex surfheat_incr(ji) = zero coldcont_incr(ji) = zero ENDDO ! ! Sum up the energy content of all layers in the soil. ! DO ji = 1, kjpindex ! IF (pcapa_en(ji,1) .LE. sn_capa) THEN ! ! Verify the energy conservation in the surface layer ! coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji)) surfheat_incr(ji) = zero ELSE ! ! Verify the energy conservation in the surface layer ! surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji)) coldcont_incr(ji) = zero ENDIF ENDDO ptn_beg(:,:) = ptn(:,:) temp_sol_beg(:) = temp_sol_new(:) END SUBROUTINE thermosoil_energy END MODULE thermosoil