| 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, inertiedat |
|---|
| 7 | use comcstfi_mod, only: pi |
|---|
| 8 | use time_phylmdz_mod, only: daysec |
|---|
| 9 | use planete_mod, only: year_day |
|---|
| 10 | use geometry_mod, only: longitude, latitude ! in radians |
|---|
| 11 | |
|---|
| 12 | implicit none |
|---|
| 13 | |
|---|
| 14 | !----------------------------------------------------------------------- |
|---|
| 15 | ! Author: Ehouarn Millour |
|---|
| 16 | ! |
|---|
| 17 | ! Purpose: Compute soil temperature using an implict 1st order scheme |
|---|
| 18 | ! |
|---|
| 19 | ! Note: depths of layers and mid-layers, soil thermal inertia and |
|---|
| 20 | ! heat capacity are commons in comsoil.h |
|---|
| 21 | !----------------------------------------------------------------------- |
|---|
| 22 | |
|---|
| 23 | !----------------------------------------------------------------------- |
|---|
| 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 | ! 0. Initialisations and preprocessing step |
|---|
| 53 | if (firstcall) then |
|---|
| 54 | ! note: firstcall is set to .true. or .false. by the caller |
|---|
| 55 | ! and not changed by soil.F |
|---|
| 56 | |
|---|
| 57 | ALLOCATE(mthermdiff(ngrid,0:nsoil-1)) ! mid-layer thermal diffusivity |
|---|
| 58 | ALLOCATE(thermdiff(ngrid,nsoil-1)) ! inter-layer thermal diffusivity |
|---|
| 59 | ALLOCATE(coefq(0:nsoil-1)) ! q_{k+1/2} coefficients |
|---|
| 60 | ALLOCATE(coefd(ngrid,nsoil-1)) ! d_k coefficients |
|---|
| 61 | ALLOCATE(alph(ngrid,nsoil-1)) ! alpha_k coefficients |
|---|
| 62 | ALLOCATE(beta(ngrid,nsoil-1)) ! beta_k coefficients |
|---|
| 63 | |
|---|
| 64 | ! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities |
|---|
| 65 | do ig=1,ngrid |
|---|
| 66 | do ik=0,nsoil-1 |
|---|
| 67 | mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa |
|---|
| 68 | !write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik) |
|---|
| 69 | enddo |
|---|
| 70 | enddo |
|---|
| 71 | |
|---|
| 72 | ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities |
|---|
| 73 | do ig=1,ngrid |
|---|
| 74 | do ik=1,nsoil-1 |
|---|
| 75 | thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik) & |
|---|
| 76 | +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1)) & |
|---|
| 77 | /(mlayer(ik)-mlayer(ik-1)) |
|---|
| 78 | ! write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik) |
|---|
| 79 | enddo |
|---|
| 80 | enddo |
|---|
| 81 | |
|---|
| 82 | ! 0.3 Build coefficients mu, q_{k+1/2}, d_k, alpha_k and capcal |
|---|
| 83 | ! mu |
|---|
| 84 | mu=mlayer(0)/(mlayer(1)-mlayer(0)) |
|---|
| 85 | |
|---|
| 86 | ! q_{1/2} |
|---|
| 87 | coefq(0)=volcapa*layer(1)/timestep |
|---|
| 88 | ! q_{k+1/2} |
|---|
| 89 | do ik=1,nsoil-1 |
|---|
| 90 | coefq(ik)=volcapa*(layer(ik+1)-layer(ik)) & |
|---|
| 91 | /timestep |
|---|
| 92 | enddo |
|---|
| 93 | |
|---|
| 94 | do ig=1,ngrid |
|---|
| 95 | ! d_k |
|---|
| 96 | do ik=1,nsoil-1 |
|---|
| 97 | coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1)) |
|---|
| 98 | enddo |
|---|
| 99 | |
|---|
| 100 | ! alph_{N-1} |
|---|
| 101 | alph(ig,nsoil-1)=coefd(ig,nsoil-1)/ & |
|---|
| 102 | (coefq(nsoil-1)+coefd(ig,nsoil-1)) |
|---|
| 103 | ! alph_k |
|---|
| 104 | do ik=nsoil-2,1,-1 |
|---|
| 105 | alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)* & |
|---|
| 106 | (1.-alph(ig,ik+1))+coefd(ig,ik)) |
|---|
| 107 | enddo |
|---|
| 108 | |
|---|
| 109 | ! capcal |
|---|
| 110 | ! Cstar |
|---|
| 111 | capcal(ig)=volcapa*layer(1)+ & |
|---|
| 112 | (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* & |
|---|
| 113 | (timestep*(1.-alph(ig,1))) |
|---|
| 114 | ! Cs |
|---|
| 115 | capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))* & |
|---|
| 116 | thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 117 | !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig) |
|---|
| 118 | enddo ! of do ig=1,ngrid |
|---|
| 119 | |
|---|
| 120 | else ! of if (firstcall) |
|---|
| 121 | |
|---|
| 122 | |
|---|
| 123 | ! 1. Compute soil temperatures |
|---|
| 124 | ! First layer: |
|---|
| 125 | do ig=1,ngrid |
|---|
| 126 | tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)* & |
|---|
| 127 | thermdiff(ig,1)/mthermdiff(ig,0))/ & |
|---|
| 128 | (1.+mu*(1.0-alph(ig,1))* & |
|---|
| 129 | thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 130 | enddo |
|---|
| 131 | ! Other layers: |
|---|
| 132 | do ik=1,nsoil-1 |
|---|
| 133 | do ig=1,ngrid |
|---|
| 134 | tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik) |
|---|
| 135 | enddo |
|---|
| 136 | enddo |
|---|
| 137 | |
|---|
| 138 | endif! of if (firstcall) |
|---|
| 139 | |
|---|
| 140 | ! 2. Compute beta coefficients (preprocessing for next time step) |
|---|
| 141 | ! Bottom layer, beta_{N-1} |
|---|
| 142 | do ig=1,ngrid |
|---|
| 143 | beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil) & |
|---|
| 144 | /(coefq(nsoil-1)+coefd(ig,nsoil-1)) |
|---|
| 145 | enddo |
|---|
| 146 | ! Other layers |
|---|
| 147 | do ik=nsoil-2,1,-1 |
|---|
| 148 | do ig=1,ngrid |
|---|
| 149 | beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+ & |
|---|
| 150 | coefd(ig,ik+1)*beta(ig,ik+1))/ & |
|---|
| 151 | (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1)) & |
|---|
| 152 | +coefd(ig,ik)) |
|---|
| 153 | enddo |
|---|
| 154 | enddo |
|---|
| 155 | |
|---|
| 156 | |
|---|
| 157 | ! 3. Compute surface diffusive flux & calorific capacity |
|---|
| 158 | do ig=1,ngrid |
|---|
| 159 | ! Cstar |
|---|
| 160 | ! capcal(ig)=volcapa(ig,1)*layer(ig,1)+ |
|---|
| 161 | ! & (thermdiff(ig,1)/(mlayer(ig,1)-mlayer(ig,0)))* |
|---|
| 162 | ! & (timestep*(1.-alph(ig,1))) |
|---|
| 163 | ! Fstar |
|---|
| 164 | |
|---|
| 165 | ! print*,'this far in soil 1' |
|---|
| 166 | ! print*,'thermdiff=',thermdiff(ig,1) |
|---|
| 167 | ! print*,'mlayer=',mlayer |
|---|
| 168 | ! print*,'beta=',beta(ig,1) |
|---|
| 169 | ! print*,'alph=',alph(ig,1) |
|---|
| 170 | ! print*,'tsoil=',tsoil(ig,1) |
|---|
| 171 | |
|---|
| 172 | fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))* & |
|---|
| 173 | (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1)) |
|---|
| 174 | |
|---|
| 175 | ! mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0)) |
|---|
| 176 | ! capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))* |
|---|
| 177 | ! & thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 178 | ! Fs |
|---|
| 179 | fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)* & |
|---|
| 180 | (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))* & |
|---|
| 181 | thermdiff(ig,1)/mthermdiff(ig,0)) & |
|---|
| 182 | -tsurf(ig)-mu*beta(ig,1)* & |
|---|
| 183 | thermdiff(ig,1)/mthermdiff(ig,0)) |
|---|
| 184 | enddo |
|---|
| 185 | |
|---|
| 186 | end |
|---|
| 187 | |
|---|
