Context Navigation

surface.F90 @ 4245

Last change on this file since 4245 was 4245, checked in by dubos, 4 years ago
simple_physics : turn zc, zd, capcal, fluxgrd into local temporaries
File size: 11.3 KB

Rev	Line
[4208]	1	MODULE surface
[4210]	2
	3	#include "use_logging.h"
	4
[4208]	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
[4228]	15	! precomputed variables
[4229]	16	REAL :: lambda
[4242]	17	REAL, ALLOCATABLE :: dz1(:),dz2(:)
[4245]	18	!$OMP THREADPRIVATE(dz1,dz2)
[4242]	19	REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
	20	!$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)
[4208]	21
[4242]	22	! internal state, written to / read from disk at checkpoint / restart
[4245]	23	REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:)
	24	!$OMP THREADPRIVATE(tsurf, tsoil)
	25	! variables below should be temporary arrays, not persistent
	26	REAL, ALLOCATABLE :: zc(:,:),zd(:,:), capcal(:), fluxgrd(:)
	27	!$OMP THREADPRIVATE(zc,zd, capcal, fluxgrd)
[4208]	28
[4245]	29	PUBLIC :: init_soil, &
	30	soil, soil_new, soil_forward, soil_backward, &
	31	zc, zd, &
[4242]	32	rnatur, albedo, emissiv, z0, inertie, &
	33	tsurf, tsoil, capcal, fluxgrd
[4228]	34
[4208]	35	CONTAINS
	36
[4245]	37	SUBROUTINE init_soil(nsoil)
	38	INTEGER, INTENT(IN) :: nsoil
[4233]	39	REAL :: min_period,dalph_soil, rk,fz1,rk1,rk2
[4208]	40	INTEGER :: jk
	41
	42	!-----------------------------------------------------------------------
[4229]	43	! ground levels
	44	! grnd=z/l where l is the skin depth of the diurnal cycle:
	45	! --------------------------------------------------------
	46
[4245]	47	WRITELOG(,) 'nsoil,firstcall=',nsoil, .TRUE.
[4210]	48
[4229]	49	ALLOCATE(dz1(nsoil),dz2(nsoil))
	50
	51	min_period=20000.
	52	dalph_soil=2.
	53
	54	! la premiere couche represente un dixieme de cycle diurne
	55	fz1=sqrt(min_period/pi)
	56
	57	DO jk=1,nsoil
	58	rk1=jk
	59	rk2=jk-1
	60	dz2(jk)=fz(rk1)-fz(rk2)
[4208]	61	ENDDO
[4229]	62	DO jk=1,nsoil-1
	63	rk1=jk+.5
	64	rk2=jk-.5
	65	dz1(jk)=1./(fz(rk1)-fz(rk2))
[4208]	66	ENDDO
[4229]	67	lambda=fz(.5)*dz1(1)
[4228]	68
[4229]	69	WRITELOG(,) 'full layers, intermediate layers (secoonds)'
	70	DO jk=1,nsoil
	71	rk=jk
	72	rk1=jk+.5
	73	rk2=jk-.5
[4210]	74	WRITELOG(,) fz(rk1)fz(rk2)pi, &
[4229]	75	& fz(rk)fz(rk)pi
[4208]	76	ENDDO
[4210]	77	LOG_INFO('init_soil')
[4228]	78
[4233]	79	CONTAINS
	80
	81	FUNCTION fz(rk) RESULT(val)
	82	REAL :: val, rk
	83	val = fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
	84	END FUNCTION fz
	85
[4208]	86	END SUBROUTINE init_soil
[4229]	87
[4245]	88	PURE SUBROUTINE soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
	89	INTEGER, INTENT(IN) :: ngrid, nsoil ! number of columns, of soil layers
	90	REAL, INTENT(IN) :: zc(ngrid, nsoil), zd(ngrid, nsoil) ! LU factorization
	91	REAL, INTENT(IN) :: ptsrf(ngrid) ! new surface temperature
	92	REAL, INTENT(INOUT) :: ptsoil(ngrid,nsoil) ! soil temperature
	93	INTEGER :: ig, jk
	94
	95	!-----------------------------------------------------------------------
	96	! Computation of the soil temperatures using the Cgrd and Dgrd
	97	! coefficient computed during the forward sweep
	98	! -----------------------------------------------
	99
	100	! surface temperature => temperature in first soil layer
	101	DO ig=1,ngrid
	102	ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ &
	103	& (lambda*(1.-zd(ig,1))+1.)
	104	ENDDO
	105
	106	! other temperatures
	107	DO jk=1,nsoil-1
	108	DO ig=1,ngrid
	109	ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
	110	ENDDO
	111	ENDDO
	112	END SUBROUTINE Soil_backward
	113
	114	PURE SUBROUTINE soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
	115	& zc, zd, pcapcal, pfluxgrd)
	116	INTEGER, INTENT(IN) :: ngrid, nsoil ! number of columns, of soil layers
	117	REAL, INTENT(IN) :: ptimestep, & ! time step
	118	& ptherm_i(ngrid), & ! thermal inertia ??
	119	& ptsrf(ngrid), & ! surface temperature before heat conduction
	120	& ptsoil(ngrid, nsoil) ! soil temperature before heat conduction
	121	REAL, INTENT(OUT) :: zc(ngrid,nsoil), &
	122	& zd(ngrid, nsoil), & ! LU factorization for backward sweep
	123	& pcapcal(ngrid), & ! effective calorific capacity
	124	& pfluxgrd(ngrid) ! conductive heat flux at the ground
	125	REAL :: z1, zdz2(ngrid)
	126	INTEGER :: jk, ig
	127
	128	!-----------------------------------------------------------------------
	129	! Computation of the Cgrd and Dgrd coefficients the backward sweep :
	130	! ---------------------------------------------------------------
	131
	132	DO jk=1,nsoil
	133	zdz2(jk)=dz2(jk)/ptimestep
	134	ENDDO
	135
	136	DO ig=1,ngrid
	137	z1=zdz2(nsoil)+dz1(nsoil-1)
	138	zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1
	139	zd(ig,nsoil-1)=dz1(nsoil-1)/z1
	140	ENDDO
	141
	142	DO jk=nsoil-1,2,-1
	143	DO ig=1,ngrid
	144	z1=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
	145	zc(ig,jk-1)= &
	146	& (ptsoil(ig,jk)zdz2(jk)+dz1(jk)zc(ig,jk))*z1
	147	zd(ig,jk-1)=dz1(jk-1)*z1
	148	ENDDO
	149	ENDDO
	150
	151	!-----------------------------------------------------------------------
	152	! computation of the surface diffusive flux from ground and
	153	! calorific capacity of the ground:
	154	! ---------------------------------
	155
	156	DO ig=1,ngrid
	157	pfluxgrd(ig)=ptherm_i(ig)dz1(1) &
	158	& (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
	159	z1=lambda*(1.-zd(ig,1))+1.
	160	pcapcal(ig)=ptherm_i(ig)* &
	161	& ptimestep(zdz2(1)+(1.-zd(ig,1))dz1(1))/z1
	162	pfluxgrd(ig)=pfluxgrd(ig) &
	163	& +pcapcal(ig)(ptsoil(ig,1)z1-lambda*zc(ig,1)-ptsrf(ig)) &
	164	& /ptimestep
	165	ENDDO
	166	END SUBROUTINE soil_forward
	167
	168	SUBROUTINE soil_new(ngrid,nsoil,ptimestep,ptherm_i, ptsrf,ptsoil, pcapcal,pfluxgrd)
	169	INTEGER, INTENT(IN) :: ngrid, nsoil ! number of columns, of soil layers
	170	REAL, INTENT(IN) :: ptimestep, & ! time step
	171	& ptherm_i(ngrid) ! thermal inertia ??
	172	REAL, INTENT(INOUT) :: ptsrf(ngrid), & ! surface temperature
	173	& ptsoil(ngrid,nsoil) ! soil temperature
	174	REAL, INTENT(OUT) :: pcapcal(ngrid), & ! effective calorific capacity
	175	& pfluxgrd(ngrid) ! conductive heat flux at the ground
	176	CALL soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
	177	CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
	178	& zc, zd, pcapcal, pfluxgrd)
	179	END SUBROUTINE soil_new
	180
[4208]	181	SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i, &
	182	& ptimestep,ptsrf,ptsoil, &
[4229]	183	& pcapcal,pfluxgrd)
	184
[4208]	185	!=======================================================================
[4229]	186	!
	187	! Auteur: Frederic Hourdin 30/01/92
	188	! -------
	189	!
	190	! objet: computation of : the soil temperature evolution
	191	! ------ the surfacic heat capacity "Capcal"
	192	! the surface conduction flux pcapcal
	193	!
	194	!
	195	! Method: implicit time integration
	196	! -------
	197	! Consecutive ground temperatures are related by:
	198	! T(k+1) = C(k) + D(k)*T(k) (1)
	199	! the coefficients C and D are computed at the t-dt time-step.
	200	! Routine structure:
	201	! 1)new temperatures are computed using (1)
	202	! 2)C and D coefficients are computed from the new temperature
	203	! profile for the t+dt time-step
	204	! 3)the coefficients A and B are computed where the diffusive
	205	! fluxes at the t+dt time-step is given by
	206	! Fdiff = A + B Ts(t+dt)
	207	! or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
	208	! with F0 = A + B (Ts(t))
	209	! Capcal = B*dt
	210	!
	211	! Interface:
	212	! ----------
	213	!
	214	! Arguments:
	215	! ----------
	216	! ngrid number of grid-points
	217	! ptimestep physical timestep (s)
	218	! pto(ngrid,nsoil) temperature at time-step t (K)
	219	! ptn(ngrid,nsoil) temperature at time step t+dt (K)
	220	! pcapcal(ngrid) specific heat (Wm-2s*K-1)
	221	! pfluxgrd(ngrid) surface diffusive flux from ground (Wm-2)
	222	!
[4208]	223	!=======================================================================
[4229]	224	! declarations:
	225	! -------------
	226
	227
[4208]	228	!-----------------------------------------------------------------------
[4229]	229	! arguments
	230	! ---------
	231
	232	INTEGER ngrid,nsoil
	233	REAL ptimestep
	234	REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid)
	235	REAL pcapcal(ngrid),pfluxgrd(ngrid)
	236	LOGICAL firstcall
	237
	238
[4208]	239	!-----------------------------------------------------------------------
[4229]	240	! local arrays
	241	! ------------
	242
	243	INTEGER ig,jk
	244	REAL zdz2(nsoil),z1(ngrid)
	245
	246	IF (firstcall) THEN
[4245]	247	! init_soil is now called by iniphyparam
	248	! CALL init_soil(ngrid, nsoil)
[4208]	249	ELSE
[4245]	250	IF(.FALSE.) THEN
	251	!-----------------------------------------------------------------------
	252	! Computation of the soil temperatures using the Cgrd and Dgrd
	253	! coefficient computed at the previous time-step:
	254	! -----------------------------------------------
	255
	256	! surface temperature
	257	DO ig=1,ngrid
	258	ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ &
	259	& (lambda*(1.-zd(ig,1))+1.)
	260	ENDDO
	261
	262	! other temperatures
	263	DO jk=1,nsoil-1
	264	DO ig=1,ngrid
	265	ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
	266	ENDDO
	267	ENDDO
	268	ELSE
	269	CALL soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
	270	END IF
	271
	272	ENDIF
	273
	274	IF(.FALSE.) THEN
[4208]	275	!-----------------------------------------------------------------------
[4245]	276	! Computation of the Cgrd and Dgrd coefficient for the next step:
	277	! ---------------------------------------------------------------
	278
	279	DO jk=1,nsoil
	280	zdz2(jk)=dz2(jk)/ptimestep
	281	ENDDO
	282
[4229]	283	DO ig=1,ngrid
[4245]	284	z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
	285	zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
	286	zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
[4208]	287	ENDDO
[4245]	288
	289	DO jk=nsoil-1,2,-1
[4229]	290	DO ig=1,ngrid
[4245]	291	z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
	292	zc(ig,jk-1)= &
	293	& (ptsoil(ig,jk)zdz2(jk)+dz1(jk)zc(ig,jk))*z1(ig)
	294	zd(ig,jk-1)=dz1(jk-1)*z1(ig)
[4208]	295	ENDDO
	296	ENDDO
[4245]	297
	298	!-----------------------------------------------------------------------
	299	! computation of the surface diffusive flux from ground and
	300	! calorific capacity of the ground:
	301	! ---------------------------------
	302
[4229]	303	DO ig=1,ngrid
[4245]	304	pfluxgrd(ig)=ptherm_i(ig)dz1(1) &
	305	& (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
	306	z1(ig)=lambda*(1.-zd(ig,1))+1.
	307	pcapcal(ig)=ptherm_i(ig)* &
	308	& ptimestep(zdz2(1)+(1.-zd(ig,1))dz1(1))/z1(ig)
	309	pfluxgrd(ig)=pfluxgrd(ig) &
	310	& +pcapcal(ig)(ptsoil(ig,1)z1(ig)-lambda*zc(ig,1)-ptsrf(ig)) &
	311	& /ptimestep
[4208]	312	ENDDO
[4245]	313	ELSE
	314	CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
	315	& zc, zd, pcapcal, pfluxgrd)
	316	END IF
[4208]	317	END SUBROUTINE soil
[4245]	318
	319	END MODULE surface
	320

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format