| 1 | MODULE surface |
|---|
| 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 | ! local saved variables: |
|---|
| 16 | ! ---------------------- |
|---|
| 17 | REAL :: lambda |
|---|
| 18 | REAL,ALLOCATABLE :: dz1(:),dz2(:),zc(:,:),zd(:,:) |
|---|
| 19 | !$OMP THREADPRIVATE(dz1,dz2,zc,zd,lambda) |
|---|
| 20 | |
|---|
| 21 | PUBLIC :: soil |
|---|
| 22 | |
|---|
| 23 | CONTAINS |
|---|
| 24 | |
|---|
| 25 | SUBROUTINE init_soil(ngrid,nsoil) |
|---|
| 26 | INTEGER, INTENT(IN) :: ngrid, nsoil |
|---|
| 27 | REAL min_period,dalph_soil |
|---|
| 28 | REAL fz,rk,fz1,rk1,rk2 |
|---|
| 29 | INTEGER :: jk |
|---|
| 30 | |
|---|
| 31 | ! this is a function definition |
|---|
| 32 | fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) |
|---|
| 33 | |
|---|
| 34 | !----------------------------------------------------------------------- |
|---|
| 35 | ! ground levels |
|---|
| 36 | ! grnd=z/l where l is the skin depth of the diurnal cycle: |
|---|
| 37 | ! -------------------------------------------------------- |
|---|
| 38 | |
|---|
| 39 | WRITELOG(*,*) 'nsoil,ngrid,firstcall=',nsoil,ngrid, .TRUE. |
|---|
| 40 | |
|---|
| 41 | ALLOCATE(dz1(nsoil),dz2(nsoil)) |
|---|
| 42 | ALLOCATE(zc(ngrid,nsoil),zd(ngrid,nsoil)) |
|---|
| 43 | |
|---|
| 44 | min_period=20000. |
|---|
| 45 | dalph_soil=2. |
|---|
| 46 | |
|---|
| 47 | ! la premiere couche represente un dixieme de cycle diurne |
|---|
| 48 | fz1=sqrt(min_period/pi) |
|---|
| 49 | |
|---|
| 50 | DO jk=1,nsoil |
|---|
| 51 | rk1=jk |
|---|
| 52 | rk2=jk-1 |
|---|
| 53 | dz2(jk)=fz(rk1)-fz(rk2) |
|---|
| 54 | ENDDO |
|---|
| 55 | DO jk=1,nsoil-1 |
|---|
| 56 | rk1=jk+.5 |
|---|
| 57 | rk2=jk-.5 |
|---|
| 58 | dz1(jk)=1./(fz(rk1)-fz(rk2)) |
|---|
| 59 | ENDDO |
|---|
| 60 | lambda=fz(.5)*dz1(1) |
|---|
| 61 | WRITELOG(*,*) 'full layers, intermediate layers (secoonds)' |
|---|
| 62 | DO jk=1,nsoil |
|---|
| 63 | rk=jk |
|---|
| 64 | rk1=jk+.5 |
|---|
| 65 | rk2=jk-.5 |
|---|
| 66 | WRITELOG(*,*) fz(rk1)*fz(rk2)*pi, & |
|---|
| 67 | & fz(rk)*fz(rk)*pi |
|---|
| 68 | ENDDO |
|---|
| 69 | |
|---|
| 70 | LOG_INFO('init_soil') |
|---|
| 71 | END SUBROUTINE init_soil |
|---|
| 72 | |
|---|
| 73 | SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i, & |
|---|
| 74 | & ptimestep,ptsrf,ptsoil, & |
|---|
| 75 | & pcapcal,pfluxgrd) |
|---|
| 76 | |
|---|
| 77 | !======================================================================= |
|---|
| 78 | ! |
|---|
| 79 | ! Auteur: Frederic Hourdin 30/01/92 |
|---|
| 80 | ! ------- |
|---|
| 81 | ! |
|---|
| 82 | ! objet: computation of : the soil temperature evolution |
|---|
| 83 | ! ------ the surfacic heat capacity "Capcal" |
|---|
| 84 | ! the surface conduction flux pcapcal |
|---|
| 85 | ! |
|---|
| 86 | ! |
|---|
| 87 | ! Method: implicit time integration |
|---|
| 88 | ! ------- |
|---|
| 89 | ! Consecutive ground temperatures are related by: |
|---|
| 90 | ! T(k+1) = C(k) + D(k)*T(k) (1) |
|---|
| 91 | ! the coefficients C and D are computed at the t-dt time-step. |
|---|
| 92 | ! Routine structure: |
|---|
| 93 | ! 1)new temperatures are computed using (1) |
|---|
| 94 | ! 2)C and D coefficients are computed from the new temperature |
|---|
| 95 | ! profile for the t+dt time-step |
|---|
| 96 | ! 3)the coefficients A and B are computed where the diffusive |
|---|
| 97 | ! fluxes at the t+dt time-step is given by |
|---|
| 98 | ! Fdiff = A + B Ts(t+dt) |
|---|
| 99 | ! or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt |
|---|
| 100 | ! with F0 = A + B (Ts(t)) |
|---|
| 101 | ! Capcal = B*dt |
|---|
| 102 | ! |
|---|
| 103 | ! Interface: |
|---|
| 104 | ! ---------- |
|---|
| 105 | ! |
|---|
| 106 | ! Arguments: |
|---|
| 107 | ! ---------- |
|---|
| 108 | ! ngrid number of grid-points |
|---|
| 109 | ! ptimestep physical timestep (s) |
|---|
| 110 | ! pto(ngrid,nsoil) temperature at time-step t (K) |
|---|
| 111 | ! ptn(ngrid,nsoil) temperature at time step t+dt (K) |
|---|
| 112 | ! pcapcal(ngrid) specific heat (W*m-2*s*K-1) |
|---|
| 113 | ! pfluxgrd(ngrid) surface diffusive flux from ground (Wm-2) |
|---|
| 114 | ! |
|---|
| 115 | !======================================================================= |
|---|
| 116 | ! declarations: |
|---|
| 117 | ! ------------- |
|---|
| 118 | |
|---|
| 119 | |
|---|
| 120 | !----------------------------------------------------------------------- |
|---|
| 121 | ! arguments |
|---|
| 122 | ! --------- |
|---|
| 123 | |
|---|
| 124 | INTEGER ngrid,nsoil |
|---|
| 125 | REAL ptimestep |
|---|
| 126 | REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid) |
|---|
| 127 | REAL pcapcal(ngrid),pfluxgrd(ngrid) |
|---|
| 128 | LOGICAL firstcall |
|---|
| 129 | |
|---|
| 130 | |
|---|
| 131 | !----------------------------------------------------------------------- |
|---|
| 132 | ! local arrays |
|---|
| 133 | ! ------------ |
|---|
| 134 | |
|---|
| 135 | INTEGER ig,jk |
|---|
| 136 | REAL za(ngrid),zb(ngrid) |
|---|
| 137 | REAL zdz2(nsoil),z1(ngrid) |
|---|
| 138 | |
|---|
| 139 | IF (firstcall) THEN |
|---|
| 140 | CALL init_soil(ngrid, nsoil) |
|---|
| 141 | ELSE |
|---|
| 142 | !----------------------------------------------------------------------- |
|---|
| 143 | ! Computation of the soil temperatures using the Cgrd and Dgrd |
|---|
| 144 | ! coefficient computed at the previous time-step: |
|---|
| 145 | ! ----------------------------------------------- |
|---|
| 146 | |
|---|
| 147 | ! surface temperature |
|---|
| 148 | DO ig=1,ngrid |
|---|
| 149 | ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ & |
|---|
| 150 | & (lambda*(1.-zd(ig,1))+1.) |
|---|
| 151 | ENDDO |
|---|
| 152 | |
|---|
| 153 | ! other temperatures |
|---|
| 154 | DO jk=1,nsoil-1 |
|---|
| 155 | DO ig=1,ngrid |
|---|
| 156 | ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) |
|---|
| 157 | ENDDO |
|---|
| 158 | ENDDO |
|---|
| 159 | |
|---|
| 160 | ENDIF |
|---|
| 161 | |
|---|
| 162 | !----------------------------------------------------------------------- |
|---|
| 163 | ! Computation of the Cgrd and Dgrd coefficient for the next step: |
|---|
| 164 | ! --------------------------------------------------------------- |
|---|
| 165 | |
|---|
| 166 | DO jk=1,nsoil |
|---|
| 167 | zdz2(jk)=dz2(jk)/ptimestep |
|---|
| 168 | ENDDO |
|---|
| 169 | |
|---|
| 170 | DO ig=1,ngrid |
|---|
| 171 | z1(ig)=zdz2(nsoil)+dz1(nsoil-1) |
|---|
| 172 | zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig) |
|---|
| 173 | zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig) |
|---|
| 174 | ENDDO |
|---|
| 175 | |
|---|
| 176 | DO jk=nsoil-1,2,-1 |
|---|
| 177 | DO ig=1,ngrid |
|---|
| 178 | z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk))) |
|---|
| 179 | zc(ig,jk-1)= & |
|---|
| 180 | & (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig) |
|---|
| 181 | zd(ig,jk-1)=dz1(jk-1)*z1(ig) |
|---|
| 182 | ENDDO |
|---|
| 183 | ENDDO |
|---|
| 184 | |
|---|
| 185 | !----------------------------------------------------------------------- |
|---|
| 186 | ! computation of the surface diffusive flux from ground and |
|---|
| 187 | ! calorific capacity of the ground: |
|---|
| 188 | ! --------------------------------- |
|---|
| 189 | |
|---|
| 190 | DO ig=1,ngrid |
|---|
| 191 | pfluxgrd(ig)=ptherm_i(ig)*dz1(1)* & |
|---|
| 192 | & (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1)) |
|---|
| 193 | z1(ig)=lambda*(1.-zd(ig,1))+1. |
|---|
| 194 | pcapcal(ig)=ptherm_i(ig)* & |
|---|
| 195 | & ptimestep*(zdz2(1)+(1.-zd(ig,1))*dz1(1))/z1(ig) |
|---|
| 196 | pfluxgrd(ig)=pfluxgrd(ig) & |
|---|
| 197 | & +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig)) & |
|---|
| 198 | & /ptimestep |
|---|
| 199 | ENDDO |
|---|
| 200 | |
|---|
| 201 | END SUBROUTINE soil |
|---|
| 202 | |
|---|
| 203 | END MODULE surface |
|---|
