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