| 1 | MODULE soil_mod |
|---|
| 2 | |
|---|
| 3 | #include "use_logging.h" |
|---|
| 4 | |
|---|
| 5 | IMPLICIT NONE |
|---|
| 6 | PRIVATE |
|---|
| 7 | SAVE |
|---|
| 8 | |
|---|
| 9 | REAL, PARAMETER :: pi=2.*ASIN(1.) |
|---|
| 10 | |
|---|
| 11 | ! common variables |
|---|
| 12 | REAL, PUBLIC :: I_mer,I_ter,Cd_mer,Cd_ter, & |
|---|
| 13 | & alb_mer,alb_ter,emi_mer,emi_ter |
|---|
| 14 | |
|---|
| 15 | ! precomputed variables |
|---|
| 16 | REAL :: lambda |
|---|
| 17 | REAL, ALLOCATABLE :: dz1(:),dz2(:) |
|---|
| 18 | !$OMP THREADPRIVATE(dz1,dz2) |
|---|
| 19 | REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:) |
|---|
| 20 | !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie) |
|---|
| 21 | |
|---|
| 22 | ! internal state, written to / read from disk at checkpoint / restart |
|---|
| 23 | REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:) |
|---|
| 24 | !$OMP THREADPRIVATE(tsurf, tsoil) |
|---|
| 25 | |
|---|
| 26 | PUBLIC :: init_soil, soil_forward, soil_backward, & |
|---|
| 27 | rnatur, albedo, emissiv, z0, inertie, & |
|---|
| 28 | tsurf, tsoil |
|---|
| 29 | |
|---|
| 30 | CONTAINS |
|---|
| 31 | |
|---|
| 32 | SUBROUTINE init_soil(nsoil) |
|---|
| 33 | INTEGER, INTENT(IN) :: nsoil |
|---|
| 34 | REAL :: min_period,dalph_soil, rk,fz1,rk1,rk2 |
|---|
| 35 | INTEGER :: jk |
|---|
| 36 | |
|---|
| 37 | !----------------------------------------------------------------------- |
|---|
| 38 | ! ground levels |
|---|
| 39 | ! grnd=z/l where l is the skin depth of the diurnal cycle: |
|---|
| 40 | ! -------------------------------------------------------- |
|---|
| 41 | |
|---|
| 42 | WRITELOG(*,*) 'nsoil,firstcall=',nsoil, .TRUE. |
|---|
| 43 | |
|---|
| 44 | ALLOCATE(dz1(nsoil),dz2(nsoil)) |
|---|
| 45 | |
|---|
| 46 | min_period=20000. |
|---|
| 47 | dalph_soil=2. |
|---|
| 48 | |
|---|
| 49 | ! la premiere couche represente un dixieme de cycle diurne |
|---|
| 50 | fz1=sqrt(min_period/pi) |
|---|
| 51 | |
|---|
| 52 | DO jk=1,nsoil |
|---|
| 53 | rk1=jk |
|---|
| 54 | rk2=jk-1 |
|---|
| 55 | dz2(jk)=fz(rk1)-fz(rk2) |
|---|
| 56 | ENDDO |
|---|
| 57 | DO jk=1,nsoil-1 |
|---|
| 58 | rk1=jk+.5 |
|---|
| 59 | rk2=jk-.5 |
|---|
| 60 | dz1(jk)=1./(fz(rk1)-fz(rk2)) |
|---|
| 61 | ENDDO |
|---|
| 62 | lambda=fz(.5)*dz1(1) |
|---|
| 63 | |
|---|
| 64 | WRITELOG(*,*) 'full layers, intermediate layers (secoonds)' |
|---|
| 65 | DO jk=1,nsoil |
|---|
| 66 | rk=jk |
|---|
| 67 | rk1=jk+.5 |
|---|
| 68 | rk2=jk-.5 |
|---|
| 69 | WRITELOG(*,*) fz(rk1)*fz(rk2)*pi, & |
|---|
| 70 | & fz(rk)*fz(rk)*pi |
|---|
| 71 | ENDDO |
|---|
| 72 | LOG_INFO('init_soil') |
|---|
| 73 | |
|---|
| 74 | CONTAINS |
|---|
| 75 | |
|---|
| 76 | FUNCTION fz(rk) RESULT(val) |
|---|
| 77 | REAL :: val, rk |
|---|
| 78 | val = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) |
|---|
| 79 | END FUNCTION fz |
|---|
| 80 | |
|---|
| 81 | END SUBROUTINE init_soil |
|---|
| 82 | |
|---|
| 83 | !======================================================================= |
|---|
| 84 | ! |
|---|
| 85 | ! Auteur: Frederic Hourdin 30/01/92 |
|---|
| 86 | ! ------- |
|---|
| 87 | ! |
|---|
| 88 | ! objet: computation of : the soil temperature evolution |
|---|
| 89 | ! ------ the surfacic heat capacity "Capcal" |
|---|
| 90 | ! the surface conduction flux pcapcal |
|---|
| 91 | ! |
|---|
| 92 | ! |
|---|
| 93 | ! Method: implicit time integration |
|---|
| 94 | ! ------- |
|---|
| 95 | ! Consecutive ground temperatures are related by: |
|---|
| 96 | ! T(k+1) = C(k) + D(k)*T(k) (1) |
|---|
| 97 | ! the coefficients C and D are computed at the t-dt time-step. |
|---|
| 98 | ! Routine structure: |
|---|
| 99 | ! 1)new temperatures are computed using (1) |
|---|
| 100 | ! 2)C and D coefficients are computed from the new temperature |
|---|
| 101 | ! profile for the t+dt time-step |
|---|
| 102 | ! 3)the coefficients A and B are computed where the diffusive |
|---|
| 103 | ! fluxes at the t+dt time-step is given by |
|---|
| 104 | ! Fdiff = A + B Ts(t+dt) |
|---|
| 105 | ! or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt |
|---|
| 106 | ! with F0 = A + B (Ts(t)) |
|---|
| 107 | ! Capcal = B*dt |
|---|
| 108 | ! |
|---|
| 109 | |
|---|
| 110 | PURE SUBROUTINE soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil) |
|---|
| 111 | INTEGER, INTENT(IN) :: ngrid, nsoil ! number of columns, of soil layers |
|---|
| 112 | REAL, INTENT(IN) :: zc(ngrid, nsoil), zd(ngrid, nsoil) ! LU factorization |
|---|
| 113 | REAL, INTENT(IN) :: ptsrf(ngrid) ! new surface temperature |
|---|
| 114 | REAL, INTENT(INOUT) :: ptsoil(ngrid,nsoil) ! soil temperature |
|---|
| 115 | INTEGER :: ig, jk |
|---|
| 116 | |
|---|
| 117 | !----------------------------------------------------------------------- |
|---|
| 118 | ! Computation of the soil temperatures using the Cgrd and Dgrd |
|---|
| 119 | ! coefficient computed during the forward sweep |
|---|
| 120 | ! ----------------------------------------------- |
|---|
| 121 | |
|---|
| 122 | ! surface temperature => temperature in first soil layer |
|---|
| 123 | DO ig=1,ngrid |
|---|
| 124 | ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ & |
|---|
| 125 | & (lambda*(1.-zd(ig,1))+1.) |
|---|
| 126 | ENDDO |
|---|
| 127 | |
|---|
| 128 | ! other temperatures |
|---|
| 129 | DO jk=1,nsoil-1 |
|---|
| 130 | DO ig=1,ngrid |
|---|
| 131 | ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) |
|---|
| 132 | ENDDO |
|---|
| 133 | ENDDO |
|---|
| 134 | END SUBROUTINE Soil_backward |
|---|
| 135 | |
|---|
| 136 | PURE SUBROUTINE soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, & |
|---|
| 137 | & zc, zd, pcapcal, pfluxgrd) |
|---|
| 138 | INTEGER, INTENT(IN) :: ngrid, nsoil ! number of columns, of soil layers |
|---|
| 139 | REAL, INTENT(IN) :: ptimestep, & ! time step |
|---|
| 140 | & ptherm_i(ngrid), & ! thermal inertia ?? |
|---|
| 141 | & ptsrf(ngrid), & ! surface temperature before heat conduction |
|---|
| 142 | & ptsoil(ngrid, nsoil) ! soil temperature before heat conduction |
|---|
| 143 | REAL, INTENT(OUT) :: zc(ngrid,nsoil), & |
|---|
| 144 | & zd(ngrid, nsoil), & ! LU factorization for backward sweep |
|---|
| 145 | & pcapcal(ngrid), & ! effective calorific capacity |
|---|
| 146 | & pfluxgrd(ngrid) ! conductive heat flux at the ground |
|---|
| 147 | REAL :: z1, zdz2(ngrid) |
|---|
| 148 | INTEGER :: jk, ig |
|---|
| 149 | |
|---|
| 150 | !----------------------------------------------------------------------- |
|---|
| 151 | ! Computation of the Cgrd and Dgrd coefficients the backward sweep : |
|---|
| 152 | ! --------------------------------------------------------------- |
|---|
| 153 | |
|---|
| 154 | DO jk=1,nsoil |
|---|
| 155 | zdz2(jk)=dz2(jk)/ptimestep |
|---|
| 156 | ENDDO |
|---|
| 157 | |
|---|
| 158 | DO ig=1,ngrid |
|---|
| 159 | z1=zdz2(nsoil)+dz1(nsoil-1) |
|---|
| 160 | zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1 |
|---|
| 161 | zd(ig,nsoil-1)=dz1(nsoil-1)/z1 |
|---|
| 162 | ENDDO |
|---|
| 163 | |
|---|
| 164 | DO jk=nsoil-1,2,-1 |
|---|
| 165 | DO ig=1,ngrid |
|---|
| 166 | z1=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk))) |
|---|
| 167 | zc(ig,jk-1)= & |
|---|
| 168 | & (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1 |
|---|
| 169 | zd(ig,jk-1)=dz1(jk-1)*z1 |
|---|
| 170 | ENDDO |
|---|
| 171 | ENDDO |
|---|
| 172 | |
|---|
| 173 | !----------------------------------------------------------------------- |
|---|
| 174 | ! computation of the surface diffusive flux from ground and |
|---|
| 175 | ! calorific capacity of the ground: |
|---|
| 176 | ! --------------------------------- |
|---|
| 177 | |
|---|
| 178 | DO ig=1,ngrid |
|---|
| 179 | pfluxgrd(ig)=ptherm_i(ig)*dz1(1)* & |
|---|
| 180 | & (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1)) |
|---|
| 181 | z1=lambda*(1.-zd(ig,1))+1. |
|---|
| 182 | pcapcal(ig)=ptherm_i(ig)* & |
|---|
| 183 | & ptimestep*(zdz2(1)+(1.-zd(ig,1))*dz1(1))/z1 |
|---|
| 184 | pfluxgrd(ig)=pfluxgrd(ig) & |
|---|
| 185 | & +pcapcal(ig)*(ptsoil(ig,1)*z1-lambda*zc(ig,1)-ptsrf(ig)) & |
|---|
| 186 | & /ptimestep |
|---|
| 187 | ENDDO |
|---|
| 188 | END SUBROUTINE soil_forward |
|---|
| 189 | |
|---|
| 190 | END MODULE soil_mod |
|---|
