| 1 | subroutine soil(ngrid,nsoil,firstcall,lastcall, |
|---|
| 2 | & therm_i, |
|---|
| 3 | & timestep,tsurf,tsoil, |
|---|
| 4 | & capcal,fluxgrd) |
|---|
| 5 | |
|---|
| 6 | use comsoil_h, only: layer, mlayer, volcapa |
|---|
| 7 | |
|---|
| 8 | implicit none |
|---|
| 9 | |
|---|
| 10 | !----------------------------------------------------------------------- |
|---|
| 11 | ! Author: Ehouarn Millour |
|---|
| 12 | ! |
|---|
| 13 | ! Purpose: Compute soil temperature using an implict 1st order scheme |
|---|
| 14 | ! |
|---|
| 15 | ! Note: depths of layers and mid-layers, soil thermal inertia and |
|---|
| 16 | ! heat capacity are commons in comsoil.h |
|---|
| 17 | !----------------------------------------------------------------------- |
|---|
| 18 | |
|---|
| 19 | ! include "dimensions.h" |
|---|
| 20 | ! include "dimphys.h" |
|---|
| 21 | |
|---|
| 22 | |
|---|
| 23 | c----------------------------------------------------------------------- |
|---|
| 24 | ! arguments |
|---|
| 25 | ! --------- |
|---|
| 26 | ! inputs: |
|---|
| 27 | integer,intent(in) :: ngrid ! number of (horizontal) grid-points |
|---|
| 28 | integer,intent(in) :: nsoil ! number of soil layers |
|---|
| 29 | logical,intent(in) :: firstcall ! identifier for initialization call |
|---|
| 30 | logical,intent(in) :: lastcall |
|---|
| 31 | real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia |
|---|
| 32 | real,intent(in) :: timestep ! time step |
|---|
| 33 | real,intent(in) :: tsurf(ngrid) ! surface temperature |
|---|
| 34 | ! outputs: |
|---|
| 35 | real,intent(out) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature |
|---|
| 36 | real,intent(out) :: capcal(ngrid) ! surface specific heat |
|---|
| 37 | real,intent(out) :: fluxgrd(ngrid) ! surface diffusive heat flux |
|---|
| 38 | |
|---|
| 39 | ! local saved variables: |
|---|
| 40 | real,dimension(:,:),save,allocatable :: mthermdiff ! mid-layer thermal diffusivity |
|---|
| 41 | real,dimension(:,:),save,allocatable :: thermdiff ! inter-layer thermal diffusivity |
|---|
| 42 | real,dimension(:),save,allocatable :: coefq ! q_{k+1/2} coefficients |
|---|
| 43 | real,dimension(:,:),save,allocatable :: coefd ! d_k coefficients |
|---|
| 44 | real,dimension(:,:),save,allocatable :: alph ! alpha_k coefficients |
|---|
| 45 | real,dimension(:,:),save,allocatable :: beta ! beta_k coefficients |
|---|
| 46 | real,save :: mu |
|---|
| 47 | !$OMP THREADPRIVATE(mthermdiff,thermdiff,coefq,coefd,alph,beta,mu) |
|---|
| 48 | |
|---|
| 49 | ! local variables: |
|---|
| 50 | integer ig,ik |
|---|
| 51 | |
|---|
| 52 | |
|---|
| 53 | ! 0. Initialisations and preprocessing step |
|---|
| 54 | if (firstcall) then |
|---|
| 55 | ! note: firstcall is set to .true. or .false. by the caller |
|---|
| 56 | ! and not changed by soil.F |
|---|
| 57 | |
|---|
| 58 | ALLOCATE(mthermdiff(ngrid,0:nsoil-1)) ! mid-layer thermal diffusivity |
|---|
| 59 | ALLOCATE(thermdiff(ngrid,nsoil-1)) ! inter-layer thermal diffusivity |
|---|
| 60 | ALLOCATE(coefq(0:nsoil-1)) ! q_{k+1/2} coefficients |
|---|
| 61 | ALLOCATE(coefd(ngrid,nsoil-1)) ! d_k coefficients |
|---|
| 62 | ALLOCATE(alph(ngrid,nsoil-1)) ! alpha_k coefficients |
|---|
| 63 | ALLOCATE(beta(ngrid,nsoil-1)) ! beta_k coefficients |
|---|
| 64 | |
|---|
| 65 | ! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities |
|---|
| 66 | do ig=1,ngrid |
|---|
| 67 | do ik=0,nsoil-1 |
|---|
| 68 | mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa |
|---|
| 69 | ! write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik) |
|---|
| 70 | enddo |
|---|
| 71 | enddo |
|---|
| 72 | |
|---|
| 73 | ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities |
|---|
| 74 | do ig=1,ngrid |
|---|
| 75 | do ik=1,nsoil-1 |
|---|
| 76 | thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik) |
|---|
| 77 | & +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1)) |
|---|
| 78 | & /(mlayer(ik)-mlayer(ik-1)) |
|---|
| 79 | ! write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik) |
|---|
| 80 | enddo |
|---|
| 81 | enddo |
|---|
| 82 | |
|---|
| 83 | ! 0.3 Build coefficients mu, q_{k+1/2}, d_k, alpha_k and capcal |
|---|
| 84 | ! mu |
|---|
| 85 | mu=mlayer(0)/(mlayer(1)-mlayer(0)) |
|---|
| 86 | |
|---|
| 87 | ! q_{1/2} |
|---|
| 88 | coefq(0)=volcapa*layer(1)/timestep |
|---|
| 89 | ! q_{k+1/2} |
|---|
| 90 | do ik=1,nsoil-1 |
|---|
| 91 | coefq(ik)=volcapa*(layer(ik+1)-layer(ik)) |
|---|
| 92 | & /timestep |
|---|
| 93 | enddo |
|---|
| 94 | |
|---|
| 95 | do ig=1,ngrid |
|---|
| 96 | ! d_k |
|---|
| 97 | do ik=1,nsoil-1 |
|---|
| 98 | coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1)) |
|---|
| 99 | enddo |
|---|
| 100 | |
|---|
| 101 | ! alph_{N-1} |
|---|
| 102 | alph(ig,nsoil-1)=coefd(ig,nsoil-1)/ |
|---|
| 103 | & (coefq(nsoil-1)+coefd(ig,nsoil-1)) |
|---|
| 104 | ! alph_k |
|---|
| 105 | do ik=nsoil-2,1,-1 |
|---|
| 106 | alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)* |
|---|
| 107 | & (1.-alph(ig,ik+1))+coefd(ig,ik)) |
|---|
| 108 | enddo |
|---|
| 109 | |
|---|
| 110 | ! capcal |
|---|
| 111 | ! Cstar |
|---|
| 112 | capcal(ig)=volcapa*layer(1)+ |
|---|
| 113 | & (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* |
|---|
| 114 | & (timestep*(1.-alph(ig,1))) |
|---|
| 115 | ! Cs |
|---|
| 116 | capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))* |
|---|
| 117 | & thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 118 | !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig) |
|---|
| 119 | enddo ! of do ig=1,ngrid |
|---|
| 120 | |
|---|
| 121 | else ! of if (firstcall) |
|---|
| 122 | |
|---|
| 123 | |
|---|
| 124 | ! 1. Compute soil temperatures |
|---|
| 125 | ! First layer: |
|---|
| 126 | do ig=1,ngrid |
|---|
| 127 | tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)* |
|---|
| 128 | & thermdiff(ig,1)/mthermdiff(ig,0))/ |
|---|
| 129 | & (1.+mu*(1.0-alph(ig,1))* |
|---|
| 130 | & thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 131 | enddo |
|---|
| 132 | ! Other layers: |
|---|
| 133 | do ik=1,nsoil-1 |
|---|
| 134 | do ig=1,ngrid |
|---|
| 135 | tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik) |
|---|
| 136 | enddo |
|---|
| 137 | enddo |
|---|
| 138 | |
|---|
| 139 | endif! of if (firstcall) |
|---|
| 140 | |
|---|
| 141 | ! 2. Compute beta coefficients (preprocessing for next time step) |
|---|
| 142 | ! Bottom layer, beta_{N-1} |
|---|
| 143 | do ig=1,ngrid |
|---|
| 144 | beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil) |
|---|
| 145 | & /(coefq(nsoil-1)+coefd(ig,nsoil-1)) |
|---|
| 146 | enddo |
|---|
| 147 | ! Other layers |
|---|
| 148 | do ik=nsoil-2,1,-1 |
|---|
| 149 | do ig=1,ngrid |
|---|
| 150 | beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+ |
|---|
| 151 | & coefd(ig,ik+1)*beta(ig,ik+1))/ |
|---|
| 152 | & (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1)) |
|---|
| 153 | & +coefd(ig,ik)) |
|---|
| 154 | enddo |
|---|
| 155 | enddo |
|---|
| 156 | |
|---|
| 157 | |
|---|
| 158 | ! 3. Compute surface diffusive flux & calorific capacity |
|---|
| 159 | do ig=1,ngrid |
|---|
| 160 | ! Cstar |
|---|
| 161 | ! capcal(ig)=volcapa(ig,1)*layer(ig,1)+ |
|---|
| 162 | ! & (thermdiff(ig,1)/(mlayer(ig,1)-mlayer(ig,0)))* |
|---|
| 163 | ! & (timestep*(1.-alph(ig,1))) |
|---|
| 164 | ! Fstar |
|---|
| 165 | |
|---|
| 166 | ! print*,'this far in soil 1' |
|---|
| 167 | ! print*,'thermdiff=',thermdiff(ig,1) |
|---|
| 168 | ! print*,'mlayer=',mlayer |
|---|
| 169 | ! print*,'beta=',beta(ig,1) |
|---|
| 170 | ! print*,'alph=',alph(ig,1) |
|---|
| 171 | ! print*,'tsoil=',tsoil(ig,1) |
|---|
| 172 | |
|---|
| 173 | fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* |
|---|
| 174 | & (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1)) |
|---|
| 175 | |
|---|
| 176 | ! mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0)) |
|---|
| 177 | ! capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))* |
|---|
| 178 | ! & thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 179 | ! Fs |
|---|
| 180 | fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)* |
|---|
| 181 | & (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))* |
|---|
| 182 | & thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 183 | & -tsurf(ig)-mu*beta(ig,1)* |
|---|
| 184 | & thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 185 | enddo |
|---|
| 186 | |
|---|
| 187 | end |
|---|
| 188 | |
|---|
