| 1 | SUBROUTINE soil(ptimestep, indice, snow, ptsrf, ptsoil, |
|---|
| 2 | s pcapcal, pfluxgrd) |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | |
|---|
| 5 | c======================================================================= |
|---|
| 6 | c |
|---|
| 7 | c Auteur: Frederic Hourdin 30/01/92 |
|---|
| 8 | c ------- |
|---|
| 9 | c |
|---|
| 10 | c objet: computation of : the soil temperature evolution |
|---|
| 11 | c ------ the surfacic heat capacity "Capcal" |
|---|
| 12 | c the surface conduction flux pcapcal |
|---|
| 13 | c |
|---|
| 14 | c |
|---|
| 15 | c Method: implicit time integration |
|---|
| 16 | c ------- |
|---|
| 17 | c Consecutive ground temperatures are related by: |
|---|
| 18 | c T(k+1) = C(k) + D(k)*T(k) (1) |
|---|
| 19 | c the coefficients C and D are computed at the t-dt time-step. |
|---|
| 20 | c Routine structure: |
|---|
| 21 | c 1)new temperatures are computed using (1) |
|---|
| 22 | c 2)C and D coefficients are computed from the new temperature |
|---|
| 23 | c profile for the t+dt time-step |
|---|
| 24 | c 3)the coefficients A and B are computed where the diffusive |
|---|
| 25 | c fluxes at the t+dt time-step is given by |
|---|
| 26 | c Fdiff = A + B Ts(t+dt) |
|---|
| 27 | c or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt |
|---|
| 28 | c with F0 = A + B (Ts(t)) |
|---|
| 29 | c Capcal = B*dt |
|---|
| 30 | c |
|---|
| 31 | c Interface: |
|---|
| 32 | c ---------- |
|---|
| 33 | c |
|---|
| 34 | c Arguments: |
|---|
| 35 | c ---------- |
|---|
| 36 | c ptimestep physical timestep (s) |
|---|
| 37 | c indice sub-surface index |
|---|
| 38 | c snow(klon,nbsrf) snow |
|---|
| 39 | c ptsrf(klon) surface temperature at time-step t (K) |
|---|
| 40 | c ptsoil(klon,nsoilmx) temperature inside the ground (K) |
|---|
| 41 | c pcapcal(klon) surfacic specific heat (W*m-2*s*K-1) |
|---|
| 42 | c pfluxgrd(klon) surface diffusive flux from ground (Wm-2) |
|---|
| 43 | c |
|---|
| 44 | c======================================================================= |
|---|
| 45 | c declarations: |
|---|
| 46 | c ------------- |
|---|
| 47 | |
|---|
| 48 | #include "dimensions.h" |
|---|
| 49 | #include "dimphy.h" |
|---|
| 50 | #include "dimsoil.h" |
|---|
| 51 | #include "indicesol.h" |
|---|
| 52 | |
|---|
| 53 | c----------------------------------------------------------------------- |
|---|
| 54 | c arguments |
|---|
| 55 | c --------- |
|---|
| 56 | |
|---|
| 57 | REAL ptimestep |
|---|
| 58 | INTEGER indice |
|---|
| 59 | REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon) |
|---|
| 60 | REAL pcapcal(klon),pfluxgrd(klon) |
|---|
| 61 | |
|---|
| 62 | c----------------------------------------------------------------------- |
|---|
| 63 | c local arrays |
|---|
| 64 | c ------------ |
|---|
| 65 | |
|---|
| 66 | INTEGER ig,jk |
|---|
| 67 | REAL zdz2(nsoilmx),z1(klon) |
|---|
| 68 | REAL min_period,dalph_soil |
|---|
| 69 | REAL ztherm_i(klon) |
|---|
| 70 | |
|---|
| 71 | c local saved variables: |
|---|
| 72 | c ---------------------- |
|---|
| 73 | REAL dz1(nsoilmx),dz2(nsoilmx) |
|---|
| 74 | REAL zc(klon,nsoilmx),zd(klon,nsoilmx) |
|---|
| 75 | REAL lambda |
|---|
| 76 | SAVE dz1,dz2,zc,zd,lambda |
|---|
| 77 | LOGICAL firstcall |
|---|
| 78 | SAVE firstcall |
|---|
| 79 | REAL isol,isno,iice |
|---|
| 80 | SAVE isol,isno,iice |
|---|
| 81 | |
|---|
| 82 | DATA firstcall/.true./ |
|---|
| 83 | |
|---|
| 84 | DATA isol,isno,iice/2000.,2000.,2000./ |
|---|
| 85 | |
|---|
| 86 | c----------------------------------------------------------------------- |
|---|
| 87 | c Depthts: |
|---|
| 88 | c -------- |
|---|
| 89 | |
|---|
| 90 | REAL fz,rk,fz1,rk1,rk2 |
|---|
| 91 | fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) |
|---|
| 92 | |
|---|
| 93 | c calcul de l'inertie thermique a partir de la variable rnat. |
|---|
| 94 | c on initialise a iice meme au-dessus d'un point de mer au cas |
|---|
| 95 | c ou le point de mer devienne point de glace au pas suivant |
|---|
| 96 | c on corrige si on a un point de terre avec ou sans glace |
|---|
| 97 | c |
|---|
| 98 | IF (indice.EQ.is_sic) THEN |
|---|
| 99 | DO ig = 1, klon |
|---|
| 100 | ztherm_i(ig) = iice |
|---|
| 101 | IF (snow(ig).GT.0.0) ztherm_i(ig) = isno |
|---|
| 102 | ENDDO |
|---|
| 103 | ELSE IF (indice.EQ.is_lic) THEN |
|---|
| 104 | DO ig = 1, klon |
|---|
| 105 | ztherm_i(ig) = iice |
|---|
| 106 | IF (snow(ig).GT.0.0) ztherm_i(ig) = isno |
|---|
| 107 | ENDDO |
|---|
| 108 | ELSE IF (indice.EQ.is_ter) THEN |
|---|
| 109 | DO ig = 1, klon |
|---|
| 110 | ztherm_i(ig) = isol |
|---|
| 111 | IF (snow(ig).GT.0.0) ztherm_i(ig) = isno |
|---|
| 112 | ENDDO |
|---|
| 113 | ELSE IF (indice.EQ.is_oce) THEN |
|---|
| 114 | DO ig = 1, klon |
|---|
| 115 | ztherm_i(ig) = iice |
|---|
| 116 | ENDDO |
|---|
| 117 | ELSE |
|---|
| 118 | PRINT*, "valeur d indice non prevue", indice |
|---|
| 119 | CALL abort |
|---|
| 120 | ENDIF |
|---|
| 121 | |
|---|
| 122 | |
|---|
| 123 | IF (firstcall) THEN |
|---|
| 124 | |
|---|
| 125 | c----------------------------------------------------------------------- |
|---|
| 126 | c ground levels |
|---|
| 127 | c grnd=z/l where l is the skin depth of the diurnal cycle: |
|---|
| 128 | c -------------------------------------------------------- |
|---|
| 129 | |
|---|
| 130 | min_period=1800. ! en secondes |
|---|
| 131 | dalph_soil=2. ! rapport entre les epaisseurs de 2 couches succ. |
|---|
| 132 | |
|---|
| 133 | OPEN(99,file='soil.def',status='old',form='formatted',err=9999) |
|---|
| 134 | READ(99,*) min_period |
|---|
| 135 | READ(99,*) dalph_soil |
|---|
| 136 | PRINT*,'Discretization for the soil model' |
|---|
| 137 | PRINT*,'First level e-folding depth',min_period, |
|---|
| 138 | s ' dalph',dalph_soil |
|---|
| 139 | CLOSE(99) |
|---|
| 140 | 9999 CONTINUE |
|---|
| 141 | |
|---|
| 142 | c la premiere couche represente un dixieme de cycle diurne |
|---|
| 143 | fz1=sqrt(min_period/3.14) |
|---|
| 144 | |
|---|
| 145 | DO jk=1,nsoilmx |
|---|
| 146 | rk1=jk |
|---|
| 147 | rk2=jk-1 |
|---|
| 148 | dz2(jk)=fz(rk1)-fz(rk2) |
|---|
| 149 | ENDDO |
|---|
| 150 | DO jk=1,nsoilmx-1 |
|---|
| 151 | rk1=jk+.5 |
|---|
| 152 | rk2=jk-.5 |
|---|
| 153 | dz1(jk)=1./(fz(rk1)-fz(rk2)) |
|---|
| 154 | ENDDO |
|---|
| 155 | lambda=fz(.5)*dz1(1) |
|---|
| 156 | PRINT*,'full layers, intermediate layers (seconds)' |
|---|
| 157 | DO jk=1,nsoilmx |
|---|
| 158 | rk=jk |
|---|
| 159 | rk1=jk+.5 |
|---|
| 160 | rk2=jk-.5 |
|---|
| 161 | PRINT *,'fz=', |
|---|
| 162 | . fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14 |
|---|
| 163 | ENDDO |
|---|
| 164 | |
|---|
| 165 | firstcall =.false. |
|---|
| 166 | |
|---|
| 167 | c Initialisations: |
|---|
| 168 | c ---------------- |
|---|
| 169 | |
|---|
| 170 | ELSE !--not firstcall |
|---|
| 171 | c----------------------------------------------------------------------- |
|---|
| 172 | c Computation of the soil temperatures using the Cgrd and Dgrd |
|---|
| 173 | c coefficient computed at the previous time-step: |
|---|
| 174 | c ----------------------------------------------- |
|---|
| 175 | |
|---|
| 176 | c surface temperature |
|---|
| 177 | DO ig=1,klon |
|---|
| 178 | ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ |
|---|
| 179 | s (lambda*(1.-zd(ig,1))+1.) |
|---|
| 180 | ENDDO |
|---|
| 181 | |
|---|
| 182 | c other temperatures |
|---|
| 183 | DO jk=1,nsoilmx-1 |
|---|
| 184 | DO ig=1,klon |
|---|
| 185 | ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) |
|---|
| 186 | ENDDO |
|---|
| 187 | ENDDO |
|---|
| 188 | |
|---|
| 189 | ENDIF !--not firstcall |
|---|
| 190 | c----------------------------------------------------------------------- |
|---|
| 191 | c Computation of the Cgrd and Dgrd coefficient for the next step: |
|---|
| 192 | c --------------------------------------------------------------- |
|---|
| 193 | |
|---|
| 194 | DO jk=1,nsoilmx |
|---|
| 195 | zdz2(jk)=dz2(jk)/ptimestep |
|---|
| 196 | ENDDO |
|---|
| 197 | |
|---|
| 198 | DO ig=1,klon |
|---|
| 199 | z1(ig)=zdz2(nsoilmx)+dz1(nsoilmx-1) |
|---|
| 200 | zc(ig,nsoilmx-1)=zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig) |
|---|
| 201 | zd(ig,nsoilmx-1)=dz1(nsoilmx-1)/z1(ig) |
|---|
| 202 | ENDDO |
|---|
| 203 | |
|---|
| 204 | DO jk=nsoilmx-1,2,-1 |
|---|
| 205 | DO ig=1,klon |
|---|
| 206 | z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk))) |
|---|
| 207 | zc(ig,jk-1)= |
|---|
| 208 | s (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig) |
|---|
| 209 | zd(ig,jk-1)=dz1(jk-1)*z1(ig) |
|---|
| 210 | ENDDO |
|---|
| 211 | ENDDO |
|---|
| 212 | |
|---|
| 213 | c----------------------------------------------------------------------- |
|---|
| 214 | c computation of the surface diffusive flux from ground and |
|---|
| 215 | c calorific capacity of the ground: |
|---|
| 216 | c --------------------------------- |
|---|
| 217 | |
|---|
| 218 | DO ig=1,klon |
|---|
| 219 | pfluxgrd(ig)=ztherm_i(ig)*dz1(1)* |
|---|
| 220 | s (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1)) |
|---|
| 221 | pcapcal(ig)=ztherm_i(ig)* |
|---|
| 222 | s (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1)) |
|---|
| 223 | z1(ig)=lambda*(1.-zd(ig,1))+1. |
|---|
| 224 | pcapcal(ig)=pcapcal(ig)/z1(ig) |
|---|
| 225 | pfluxgrd(ig)=pfluxgrd(ig) |
|---|
| 226 | s +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig)) |
|---|
| 227 | s /ptimestep |
|---|
| 228 | ENDDO |
|---|
| 229 | |
|---|
| 230 | RETURN |
|---|
| 231 | END |
|---|
