--- trunk/Sources/phylmd/Interface_surf/soil.f 2015/07/20 16:01:49 157 +++ trunk/Sources/phylmd/Interface_surf/soil.f 2016/02/05 16:02:34 175 @@ -4,76 +4,68 @@ contains - SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil, pcapcal, & - pfluxgrd) + SUBROUTINE soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux) ! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1 2004/05/19 - USE dimens_m - USE indicesol - USE dimphy - USE dimsoil - USE suphec_m - - ! ======================================================================= - - ! Auteur: Frederic Hourdin 30/01/92 - ! ------- - - ! objet: computation of : the soil temperature evolution - ! ------ the surfacic heat capacity "Capcal" - ! the surface conduction flux pcapcal + ! Author: Frederic Hourdin 30/01/92 + ! Object: computation of the soil temperature evolution, the + ! surfacic heat capacity "Soilcap" and the surface conduction flux ! Method: implicit time integration - ! ------- + ! Consecutive ground temperatures are related by: - ! T(k+1) = C(k) + D(k)*T(k) (1) + ! T(k+1) = C(k) + D(k)*T(k) (1) ! the coefficients C and D are computed at the t-dt time-step. ! Routine structure: - ! 1)new temperatures are computed using (1) - ! 2)C and D coefficients are computed from the new temperature + ! 1) new temperatures are computed using (1) + ! 2) C and D coefficients are computed from the new temperature ! profile for the t+dt time-step - ! 3)the coefficients A and B are computed where the diffusive + ! 3) the coefficients A and B are computed where the diffusive ! fluxes at the t+dt time-step is given by ! Fdiff = A + B Ts(t+dt) - ! or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt + ! or Fdiff = F0 + Soilcap (Ts(t+dt)-Ts(t))/dt ! with F0 = A + B (Ts(t)) - ! Capcal = B*dt + ! Soilcap = B*dt + + USE dimens_m, only: + USE indicesol + USE dimphy + USE dimsoil + USE suphec_m ! Interface: ! ---------- ! Arguments: ! ---------- - ! ptimestep physical timestep (s) - ! indice sub-surface index - ! snow(klon,nbsrf) snow - ! ptsrf(knon) surface temperature at time-step t (K) - ! ptsoil(klon,nsoilmx) temperature inside the ground (K) - ! pcapcal(klon) surfacic specific heat (W*m-2*s*K-1) - ! pfluxgrd(klon) surface diffusive flux from ground (Wm-2) + ! dtime physical timestep (s) + ! indice sub-surface index + ! snow(klon, nbsrf) snow + ! tsurf(knon) surface temperature at time-step t (K) + ! tsoil(klon, nsoilmx) temperature inside the ground (K) + ! soilcap(klon) surfacic specific heat (W*m-2*s*K-1) + ! soilflux(klon) surface diffusive flux from ground (Wm-2) - ! ======================================================================= ! declarations: ! ------------- - ! ----------------------------------------------------------------------- ! arguments ! --------- - REAL ptimestep - INTEGER indice, knon - REAL ptsrf(knon), ptsoil(klon, nsoilmx), snow(klon) - REAL pcapcal(klon), pfluxgrd(klon) + REAL dtime + INTEGER nisurf, knon + REAL tsurf(knon), tsoil(klon, nsoilmx), snow(klon) + REAL soilcap(klon), soilflux(klon) ! ----------------------------------------------------------------------- ! local arrays ! ------------ INTEGER ig, jk - ! $$$ REAL zdz2(nsoilmx),z1(klon) + ! $$$ REAL zdz2(nsoilmx), z1(klon) REAL zdz2(nsoilmx), z1(klon, nbsrf) REAL min_period, dalph_soil REAL ztherm_i(klon) @@ -81,7 +73,7 @@ ! local saved variables: ! ---------------------- REAL dz1(nsoilmx), dz2(nsoilmx) - ! $$$ REAL zc(klon,nsoilmx),zd(klon,nsoilmx) + ! $$$ REAL zc(klon, nsoilmx), zd(klon, nsoilmx) REAL zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf) REAL lambda SAVE dz1, dz2, zc, zd, lambda @@ -100,37 +92,37 @@ REAL rk, fz1, rk1, rk2 - pfluxgrd(:) = 0. + soilflux(:) = 0. ! calcul de l'inertie thermique a partir de la variable rnat. ! on initialise a iice meme au-dessus d'un point de mer au cas ! ou le point de mer devienne point de glace au pas suivant ! on corrige si on a un point de terre avec ou sans glace - IF (indice==is_sic) THEN + IF (nisurf==is_sic) THEN DO ig = 1, knon ztherm_i(ig) = iice IF (snow(ig)>0.0) ztherm_i(ig) = isno END DO - ELSE IF (indice==is_lic) THEN + ELSE IF (nisurf==is_lic) THEN DO ig = 1, knon ztherm_i(ig) = iice IF (snow(ig)>0.0) ztherm_i(ig) = isno END DO - ELSE IF (indice==is_ter) THEN + ELSE IF (nisurf==is_ter) THEN DO ig = 1, knon ztherm_i(ig) = isol IF (snow(ig)>0.0) ztherm_i(ig) = isno END DO - ELSE IF (indice==is_oce) THEN + ELSE IF (nisurf==is_oce) THEN DO ig = 1, knon ztherm_i(ig) = iice END DO ELSE - PRINT *, 'valeur d indice non prevue', indice + PRINT *, 'valeur d indice non prevue', nisurf STOP 1 END IF - IF (firstsurf(indice)) THEN + IF (firstsurf(nisurf)) THEN ! ----------------------------------------------------------------------- ! ground levels @@ -144,7 +136,7 @@ READ (99, *) min_period READ (99, *) dalph_soil PRINT *, 'Discretization for the soil model' - PRINT *, 'First level e-folding depth', min_period, ' dalph', & + PRINT *, 'First level e-folding depth', min_period, ' dalph', & dalph_soil CLOSE (99) 9999 CONTINUE @@ -171,7 +163,7 @@ PRINT *, 'fz=', fz(rk1)*fz(rk2)*3.14, fz(rk)*fz(rk)*3.14 END DO ! PB - firstsurf(indice) = .FALSE. + firstsurf(nisurf) = .FALSE. ! Initialisations: ! ---------------- @@ -184,14 +176,14 @@ ! surface temperature DO ig = 1, knon - ptsoil(ig, 1) = (lambda*zc(ig,1,indice)+ptsrf(ig))/(lambda*(1.-zd(ig,1, & - indice))+1.) + tsoil(ig, 1) = (lambda*zc(ig, 1, nisurf)+tsurf(ig))/(lambda*(1.-zd(ig, 1, & + nisurf))+1.) END DO ! other temperatures DO jk = 1, nsoilmx - 1 DO ig = 1, knon - ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, & + tsoil(ig, jk+1) = zc(ig, jk, nisurf) + zd(ig, jk, nisurf)*tsoil(ig, & jk) END DO END DO @@ -201,30 +193,30 @@ ! Computation of the Cgrd and Dgrd coefficient for the next step: ! --------------------------------------------------------------- - ! $$$ PB ajout pour cas glace de mer - IF (indice==is_sic) THEN + ! $$$ PB ajout pour cas glace de mer + IF (nisurf==is_sic) THEN DO ig = 1, knon - ptsoil(ig, nsoilmx) = rtt - 1.8 + tsoil(ig, nsoilmx) = rtt - 1.8 END DO END IF DO jk = 1, nsoilmx - zdz2(jk) = dz2(jk)/ptimestep + zdz2(jk) = dz2(jk)/dtime END DO DO ig = 1, knon - z1(ig, indice) = zdz2(nsoilmx) + dz1(nsoilmx-1) - zc(ig, nsoilmx-1, indice) = zdz2(nsoilmx)*ptsoil(ig, nsoilmx)/ & - z1(ig, indice) - zd(ig, nsoilmx-1, indice) = dz1(nsoilmx-1)/z1(ig, indice) + z1(ig, nisurf) = zdz2(nsoilmx) + dz1(nsoilmx-1) + zc(ig, nsoilmx-1, nisurf) = zdz2(nsoilmx)*tsoil(ig, nsoilmx)/ & + z1(ig, nisurf) + zd(ig, nsoilmx-1, nisurf) = dz1(nsoilmx-1)/z1(ig, nisurf) END DO DO jk = nsoilmx - 1, 2, -1 DO ig = 1, knon - z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice))) - zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) & - )*z1(ig, indice) - zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice) + z1(ig, nisurf) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig, jk, nisurf))) + zc(ig, jk-1, nisurf) = (tsoil(ig, jk)*zdz2(jk)+dz1(jk)*zc(ig, jk, nisurf) & + )*z1(ig, nisurf) + zd(ig, jk-1, nisurf) = dz1(jk-1)*z1(ig, nisurf) END DO END DO @@ -234,13 +226,13 @@ ! --------------------------------- DO ig = 1, knon - pfluxgrd(ig) = ztherm_i(ig)*dz1(1)*(zc(ig,1,indice)+(zd(ig,1, & - indice)-1.)*ptsoil(ig,1)) - pcapcal(ig) = ztherm_i(ig)*(dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1)) - z1(ig, indice) = lambda*(1.-zd(ig,1,indice)) + 1. - pcapcal(ig) = pcapcal(ig)/z1(ig, indice) - pfluxgrd(ig) = pfluxgrd(ig) + pcapcal(ig)*(ptsoil(ig,1)*z1(ig,indice)- & - lambda*zc(ig,1,indice)-ptsrf(ig))/ptimestep + soilflux(ig) = ztherm_i(ig)*dz1(1)*(zc(ig, 1, nisurf)+(zd(ig, 1, & + nisurf)-1.)*tsoil(ig, 1)) + soilcap(ig) = ztherm_i(ig)*(dz2(1)+dtime*(1.-zd(ig, 1, nisurf))*dz1(1)) + z1(ig, nisurf) = lambda*(1.-zd(ig, 1, nisurf)) + 1. + soilcap(ig) = soilcap(ig)/z1(ig, nisurf) + soilflux(ig) = soilflux(ig) + soilcap(ig)*(tsoil(ig, 1)*z1(ig, nisurf)- & + lambda*zc(ig, 1, nisurf)-tsurf(ig))/dtime END DO contains