Changeset 4245 for dynamico_lmdz
- Timestamp:
- Jun 24, 2020, 11:27:16 AM (4 years ago)
- Location:
- dynamico_lmdz/simple_physics
- Files:
-
- 7 edited
-
config/DYNAMICO/TEST/file_def_physics.xml (modified) (2 diffs)
-
config/DYNAMICO/TEST/job_CICLAD.sh (modified) (1 diff)
-
config/DYNAMICO/TEST/run.def (modified) (1 diff)
-
phyparam/physics/iniphyparam_mod.F90 (modified) (3 diffs)
-
phyparam/physics/phyparam_mod.F90 (modified) (12 diffs)
-
phyparam/physics/radiative_sw.F90 (modified) (1 diff)
-
phyparam/physics/surface.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/file_def_physics.xml
r4244 r4245 10 10 <field field_ref="phyparam_swtop" /> 11 11 <field field_ref="phyparam_lwsurf" /> 12 <field field_ref="ps"/> 12 13 13 <field field_ref="phyparam_u" />14 <field field_ref="phyparam_v" />15 14 <field field_ref="phyparam_temp" /> 16 <field field_ref="phyparam_du" />17 <field field_ref="phyparam_dv" />18 <field field_ref="phyparam_dtsw" />19 <field field_ref="phyparam_dtlw" />20 15 <field field_ref="phyparam_swflux_down" /> 21 16 <field field_ref="phyparam_swflux_up" /> … … 23 18 <field field_ref="phyparam_lwflux_up" /> 24 19 25 </field_group> 20 <!-- 21 <field field_ref="phyparam_u" /> 22 <field field_ref="phyparam_v" /> 23 <field field_ref="phyparam_du" /> 24 <field field_ref="phyparam_dv" /> 25 <field field_ref="phyparam_dtsw" /> 26 <field field_ref="phyparam_dtlw" /> 27 --> 28 29 </field_group> 26 30 27 31 <variable name="model" type="string" > dynamico </variable> -
dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/job_CICLAD.sh
r4239 r4245 4 4 #PBS -n 5 5 #PBS -l nodes=1:ppn=60 6 #PBS -l walltime=0 1:00:006 #PBS -l walltime=04:00:00 7 7 #PBS -l mem=120gb 8 8 #PBS -l vmem=120gb -
dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/run.def
r4242 r4245 42 42 43 43 #---------------- Run --------------- 44 run_length=2592000 44 run_length=8640000 45 #run_length=2592000 45 46 # 21600s = 6h 46 47 write_period=21600 -
dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90
r4241 r4245 71 71 72 72 SUBROUTINE iniphyparam(ptimestep, punjours, prad, pg, pr, pcpp) 73 USE comgeomfi, ONLY : nsoilmx 74 USE surface, ONLY : init_soil 73 75 USE phys_const, ONLY : planet_rad,g,r,cpp,rcp,dtphys,unjours 74 76 USE callkeys … … 76 78 77 79 CALL read_params(ptimestep) 78 79 80 CALL check_mismatch('unjours', punjours, unjours) 80 81 CALL check_mismatch('rad', prad, planet_rad) … … 98 99 99 100 LOG_INFO('iniphyparam') 101 102 CALL init_soil(nsoilmx) 100 103 END SUBROUTINE iniphyparam 101 104 -
dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90
r4242 r4245 3 3 USE callkeys 4 4 USE comgeomfi 5 USE surface6 5 IMPLICIT NONE 7 6 PRIVATE … … 25 24 & pdu,pdv,pdt,pdpsrf) BIND(C, name='phyparam_phyparam') 26 25 USE phys_const, ONLY : g, rcp, r, unjours 26 USE surface, ONLY : soil_forward, soil_backward 27 USE surface, ONLY : z0, inertie, emissiv, albedo ! precomputed 28 USE surface, ONLY : tsurf, tsoil ! state variables 29 ! USE surface, ONLY : capcal, fluxgrd, zc, zd ! should be temporaries 27 30 USE turbulence, ONLY : vdif 28 31 USE convection, ONLY : convadj … … 66 69 67 70 ! Local variables : 68 REAL, DIMENSION(ngrid) :: mu0 69 INTEGER :: j,l,ig,nlevel,igout 70 LOGICAL :: lwrite 71 ! 72 REAL :: zday, zdtime 73 REAL zh(ngrid,nlayer),z1,z2 74 REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer) 75 REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer) 76 REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid) 77 REAL zflubid(ngrid),zpmer(ngrid) 78 REAL zpopsk(ngrid,nlayer) 71 REAL :: zh(ngrid,nlayer), & ! potential temperature 72 & zpopsk(ngrid,nlayer), & ! Exner function 73 & zzlev(ngrid,nlayer+1), & ! altitude of interfaces 74 & zzlay(ngrid,nlayer), & ! altitude of full levels 75 & fluxrad(ngrid), & ! radiative flux at surface 76 & zc(ngrid, nsoilmx), & ! LU coefficients for soil implicit solve 77 & zd(ngrid, nsoilmx), & 78 & fluxgrd(ngrid), & ! heat flux from deep soil 79 & capcal(ngrid) & ! effective heat capacity of soil 80 & zdufr(ngrid,nlayer), & ! partial tendencies for zonal velocity, 81 & zdvfr(ngrid,nlayer), & ! meridional velocity, 82 & zdhfr(ngrid,nlayer), & ! potential temperature, 83 & zdtsrfr(ngrid), & ! surface temperature 84 & zdtsrf(ngrid), & ! total tendency of surface temperature 85 & zflubid(ngrid), & ! radiative + deep soil fluxes 86 & zpmer(ngrid), & ! sea-level pressure 79 87 REAL zdum1(ngrid,nlayer) 80 88 REAL zdum2(ngrid,nlayer) 81 89 REAL zdum3(ngrid,nlayer) 82 REAL fluxrad(ngrid) 90 91 INTEGER :: j,l,ig,nlevel,igout 92 LOGICAL :: lwrite 93 REAL :: zday, zdtime, z1, z2 83 94 84 95 WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid) … … 88 99 IF (ngrid.NE.ngridmax) THEN 89 100 PRINT*,'STOP in inifis' 90 PRINT*,'Probleme de dimen esions :'101 PRINT*,'Probleme de dimensions :' 91 102 PRINT*,'ngrid = ',ngrid 92 103 PRINT*,'ngridmax = ',ngridmax … … 159 170 ENDDO 160 171 172 !------------------------------------------------------------- 173 ! soil temperatures : 1st half of implicit time integration 174 ! forward sweep from deep ground to surface 175 ! yields LU coefficients zc,zd and capcal, fluxgrd 176 ! ---------------------------------------------------------- 177 178 IF (callsoil) THEN 179 CALL soil_forward(ngrid,nsoilmx, ptimestep, inertie, tsurf, tsoil, & 180 & zc, zd, capcal, fluxgrd) 181 182 ! CALL soil_new(ngrid,nsoilmx,ptimestep,inertie, & 183 ! tsurf, tsoil, capcal,fluxgrd) 184 ! CALL soil(ngrid,nsoilmx,.false.,inertie, & 185 ! & ptimestep,tsurf,tsoil,capcal,fluxgrd) 186 ELSE 187 capcal(:) = capcal_nosoil 188 fluxgrd(:) = 0. 189 ENDIF 190 191 IF(lverbose) THEN 192 WRITELOG(*,*) 'Surface Heat capacity, conduction Flux, Ts' 193 WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout) 194 LOG_DBG('phyparam') 195 ENDIF 196 161 197 !----------------------------------------------------------------------- 162 198 ! 2. Compute radiative tendencies : … … 169 205 !----------------------------------------------------------------------- 170 206 ! 3. Vertical diffusion (turbulent mixing): 171 ! ----------------------------------------- 207 ! Kz is computed then vertical diffusion is integrated in time implicitly 208 ! using a linear relationship between surface heat flux and air temperature 209 ! in lowest level (Robin-type BC) 210 ! ------------------------------------------------------------------- 172 211 ! 173 212 IF(calldifv) THEN … … 212 251 ENDDO 213 252 ENDIF 253 254 !------------------------------------------------------------- 255 ! soil temperatures : 2nd half of implicit time integration 256 ! using updated tsurf as input 257 ! ---------------------------------------------------------- 258 259 DO ig=1,ngrid 260 tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) 261 ENDDO 262 263 WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1) 264 265 IF (callsoil) THEN 266 CALL soil_backward(ngrid,nsoilmx, zc,zd, tsurf,tsoil) 267 IF(lverbose) THEN 268 WRITELOG(*,*) 'Surface Ts, dTs, dt' 269 WRITELOG(*,*) tsurf(igout), zdtsrf(igout), ptimestep 270 LOG_DBG('phyparam') 271 ENDIF 272 END IF 273 274 214 275 ! 215 276 !----------------------------------------------------------------------- … … 243 304 ENDDO 244 305 245 ENDIF246 247 !-----------------------------------------------------------------------248 ! On ajoute les tendances physiques a la temperature du sol:249 ! ---------------------------------------------------------------250 251 DO ig=1,ngrid252 tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)253 ENDDO254 255 WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)256 257 !-----------------------------------------------------------------------258 ! soil temperatures:259 ! --------------------260 261 IF (callsoil) THEN262 CALL soil(ngrid,nsoilmx,.false.,inertie, &263 & ptimestep,tsurf,tsoil,capcal,fluxgrd)264 IF(lverbose) THEN265 WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt'266 WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout), &267 & zdtsrf(igout), ptimestep268 LOG_DBG('phyparam')269 ENDIF270 306 ENDIF 271 307 … … 309 345 !$cython wrapper def alloc(ngrid, nlayer) : phy.phyparam_alloc(ngrid, nlayer) 310 346 USE astronomy, ONLY : iniorbit 347 USE surface 311 348 INTEGER, INTENT(IN), VALUE :: ngrid, nlayer 349 ! allocate precomputed arrays 350 ALLOCATE(rnatur(ngrid), albedo(ngrid), emissiv(ngrid)) 351 ALLOCATE(z0(ngrid),inertie(ngrid)) 312 352 ! allocate arrays for internal state 313 353 ALLOCATE(tsurf(ngrid)) 314 354 ALLOCATE(tsoil(ngrid,nsoilmx)) 315 ! we could avoid the arrays below with a different implementation of surface / radiation / turbulence coupling 316 ALLOCATE(capcal(ngrid),fluxgrd(ngrid)) 317 ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx)) 318 ! allocate precomputed arrays 319 ALLOCATE(rnatur(ngrid), albedo(ngrid), emissiv(ngrid)) 320 ALLOCATE(z0(ngrid),inertie(ngrid)) 355 IF(.FALSE.) THEN ! arrays below are now local temporaries in phyparam 356 ALLOCATE(capcal(ngrid),fluxgrd(ngrid)) 357 ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx)) 358 END IF 321 359 CALL iniorbit 322 360 END SUBROUTINE alloc … … 326 364 !$cython wrapper def precompute() : phy.phyparam_precompute() 327 365 ! precompute time-independent arrays 366 USE surface 328 367 rnatur(:) = 1. 329 368 inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter … … 337 376 !$cython wrapper def coldstart (ngrid, timestep): phy.phyparam_coldstart(ngrid, timestep) 338 377 ! create internal state to start a run without a restart file 378 USE surface 339 379 INTEGER, INTENT(IN), VALUE :: ngrid 340 380 REAL, INTENT(IN), VALUE :: ptimestep … … 343 383 icount=0 344 384 IF(callsoil) THEN 345 ! initializes zc, zd, capcal, fluxgrd 346 CALL soil(ngrid,nsoilmx,.TRUE.,inertie, & 347 & ptimestep,tsurf,tsoil,capcal,fluxgrd) 385 IF(.FALSE.) THEN 386 ! init_soil is now called by iniphyparam 387 ! initializes zc, zd, capcal, fluxgrd 388 CALL soil(ngrid,nsoilmx,.TRUE.,inertie, & 389 & ptimestep,tsurf,tsoil,capcal,fluxgrd) 390 END IF 391 IF(.FALSE.) THEN ! soil_forward is now called by phyparam 392 CALL soil_forward(ngrid, nsoilmx, ptimestep, inertie, tsurf, tsoil, & 393 & zc, zd, capcal, fluxgrd) 394 END IF 348 395 ELSE 349 396 WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off' -
dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90
r4243 r4245 69 69 REAL :: buf1(ngrid), buf2(ngrid, nlayer+1) ! buffers for output 70 70 ! fluxes are non-zero only on those points where the sun shines (mu0>0) 71 ! We compute only on those ncount points and gather them to vectorize 71 ! We compute only on those ncount points and gather them to vectorize 72 72 INTEGER :: ncount, index(ngrid) 73 73 ! In the work arrays below, ngrid should be ncount but ncount is not known yet -
dynamico_lmdz/simple_physics/phyparam/physics/surface.F90
r4242 r4245 16 16 REAL :: lambda 17 17 REAL, ALLOCATABLE :: dz1(:),dz2(:) 18 !$OMP THREADPRIVATE(dz1,dz2 ,zc)18 !$OMP THREADPRIVATE(dz1,dz2) 19 19 REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:) 20 20 !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie) 21 21 22 22 ! internal state, written to / read from disk at checkpoint / restart 23 REAL, ALLOCATABLE :: zc(:,:),zd(:,:) 24 !$OMP THREADPRIVATE(zc,zd) 25 REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:), capcal(:), fluxgrd(:) 26 !$OMP THREADPRIVATE(tsurf, tsoil, capcal, fluxgrd) 27 28 PUBLIC :: soil, zc, zd, & 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) 28 29 PUBLIC :: init_soil, & 30 soil, soil_new, soil_forward, soil_backward, & 31 zc, zd, & 29 32 rnatur, albedo, emissiv, z0, inertie, & 30 33 tsurf, tsoil, capcal, fluxgrd … … 32 35 CONTAINS 33 36 34 SUBROUTINE init_soil(n grid,nsoil)35 INTEGER, INTENT(IN) :: n grid, nsoil37 SUBROUTINE init_soil(nsoil) 38 INTEGER, INTENT(IN) :: nsoil 36 39 REAL :: min_period,dalph_soil, rk,fz1,rk1,rk2 37 40 INTEGER :: jk … … 42 45 ! -------------------------------------------------------- 43 46 44 WRITELOG(*,*) 'nsoil, ngrid,firstcall=',nsoil,ngrid, .TRUE.47 WRITELOG(*,*) 'nsoil,firstcall=',nsoil, .TRUE. 45 48 46 49 ALLOCATE(dz1(nsoil),dz2(nsoil)) … … 82 85 83 86 END SUBROUTINE init_soil 87 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 84 180 85 181 SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i, & … … 149 245 150 246 IF (firstcall) THEN 151 CALL init_soil(ngrid, nsoil) 247 ! init_soil is now called by iniphyparam 248 ! CALL init_soil(ngrid, nsoil) 152 249 ELSE 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 153 275 !----------------------------------------------------------------------- 154 ! Computation of the soil temperatures using the Cgrd and Dgrd 155 ! coefficient computed at the previous time-step: 156 ! ----------------------------------------------- 157 158 ! surface temperature 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 159 283 DO ig=1,ngrid 160 ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ &161 & (lambda*(1.-zd(ig,1))+1.)162 ENDDO163 164 ! other temperatures165 DO jk= 1,nsoil-1284 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) 287 ENDDO 288 289 DO jk=nsoil-1,2,-1 166 290 DO ig=1,ngrid 167 ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) 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) 168 295 ENDDO 169 296 ENDDO 170 171 ENDIF 172 173 !----------------------------------------------------------------------- 174 ! Computation of the Cgrd and Dgrd coefficient for the next step: 175 ! --------------------------------------------------------------- 176 177 DO jk=1,nsoil 178 zdz2(jk)=dz2(jk)/ptimestep 179 ENDDO 180 181 DO ig=1,ngrid 182 z1(ig)=zdz2(nsoil)+dz1(nsoil-1) 183 zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig) 184 zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig) 185 ENDDO 186 187 DO jk=nsoil-1,2,-1 297 298 !----------------------------------------------------------------------- 299 ! computation of the surface diffusive flux from ground and 300 ! calorific capacity of the ground: 301 ! --------------------------------- 302 188 303 DO ig=1,ngrid 189 z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk))) 190 zc(ig,jk-1)= & 191 & (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig) 192 zd(ig,jk-1)=dz1(jk-1)*z1(ig) 193 ENDDO 194 ENDDO 195 196 !----------------------------------------------------------------------- 197 ! computation of the surface diffusive flux from ground and 198 ! calorific capacity of the ground: 199 ! --------------------------------- 200 201 DO ig=1,ngrid 202 pfluxgrd(ig)=ptherm_i(ig)*dz1(1)* & 203 & (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1)) 204 z1(ig)=lambda*(1.-zd(ig,1))+1. 205 pcapcal(ig)=ptherm_i(ig)* & 206 & ptimestep*(zdz2(1)+(1.-zd(ig,1))*dz1(1))/z1(ig) 207 pfluxgrd(ig)=pfluxgrd(ig) & 208 & +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig)) & 209 & /ptimestep 210 ENDDO 211 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 312 ENDDO 313 ELSE 314 CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, & 315 & zc, zd, pcapcal, pfluxgrd) 316 END IF 212 317 END SUBROUTINE soil 213 214 END MODULE surface 318 319 END MODULE surface 320
Note: See TracChangeset
for help on using the changeset viewer.
