Changeset 4228 for dynamico_lmdz/simple_physics
- Timestamp:
- Jan 15, 2020, 11:43:11 PM (5 years ago)
- Location:
- dynamico_lmdz/simple_physics/phyparam
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/phyparam/Makefile
r4227 r4228 30 30 obj/convection.o : logging.o 31 31 obj/surface.o : logging.o 32 obj/astronomy.o : logging.o planet.o32 obj/astronomy.o : logging.o 33 33 obj/solar.o : logging.o astronomy.o 34 34 obj/radiative_lw.o : logging.o phys_const.o 35 35 obj/radiative_sw.o : logging.o phys_const.o 36 obj/turbulence.o : logging.o phys_const.o planet.o37 MAIN : callkeys.o phys_const.o astronomy.o turbulence.o surface.o convection.o 36 obj/turbulence.o : logging.o phys_const.o 37 MAIN : callkeys.o phys_const.o astronomy.o turbulence.o surface.o convection.o planet.o 38 38 obj/iniphyparam_mod.o : MAIN read_param_mod.o 39 39 obj/phyparam_mod.o : MAIN solar.o radiative_sw.o radiative_lw.o comgeomfi.o writefield_mod.o -
dynamico_lmdz/simple_physics/phyparam/dynphy_lonlat/physiq_mod.F90
r4223 r4228 82 82 d_u,d_v,d_t,d_ps) 83 83 84 d_qx(:,:,:)=0. 85 84 86 IF(lafin) THEN 85 87 call iotd_fin -
dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90
r4210 r4228 14 14 15 15 SUBROUTINE solarlong(pday,psollong) 16 USE planet17 16 REAL, INTENT(IN) :: pday ! jour de l'annee (le jour 0 correspondant a l'equinoxe) 18 17 REAL, INTENT(OUT) :: psollong ! solar longitude -
dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90
r4223 r4228 1 1 MODULE iniphyparam_mod 2 #include "use_logging.h" 2 3 IMPLICIT NONE 3 4 PRIVATE … … 27 28 REAL, INTENT(IN) :: ptimestep, prad, pg, pr, pcpp, punjours 28 29 29 print*,'INIPHYPARAM'30 print*,'Avant les getpar '31 32 30 CALL read_param('unjours', 86400., unjours,'unjours') 33 31 CALL read_param('planet_rad',prad,planet_rad,'planet_rad') … … 72 70 CALL check_mismatch('cpp', pcpp, cpp) 73 71 74 print*,'Activation de la physique:'75 print*,' R=',r76 print*,' Rayonnement ',callrad77 print*,' Diffusion verticale turbulente ', calldifv78 print*,' Ajustement convectif ',calladj79 print*,' Sol ',callsoil80 print*,' Cycle diurne ',diurnal81 72 WRITELOG(*,*) 'Activation de la physique:' 73 WRITELOG(*,*) ' R=',r 74 WRITELOG(*,*) ' Rayonnement ',callrad 75 WRITELOG(*,*) ' Diffusion verticale turbulente ', calldifv 76 WRITELOG(*,*) ' Ajustement convectif ',calladj 77 WRITELOG(*,*) ' Sol ',callsoil 78 WRITELOG(*,*) ' Cycle diurne ',diurnal 79 82 80 ! choice of the frequency of the computation of radiations 83 81 IF(diurnal) THEN … … 87 85 ENDIF 88 86 iradia=1 89 PRINT*,'unjours',punjours90 PRINT*,'The radiative transfer is computed each ', &87 WRITELOG(*,*) 'unjours',punjours 88 WRITELOG(*,*) 'The radiative transfer is computed each ', & 91 89 & iradia,' physical time-step or each ', & 92 90 & iradia*ptimestep,' seconds' … … 94 92 dtphys=ptimestep 95 93 94 LOG_INFO('iniphyparama') 96 95 END SUBROUTINE iniphyparam 97 96 -
dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90
r4223 r4228 1 1 MODULE phyparam_mod 2 #include "use_logging.h" 2 3 USE callkeys 3 4 USE comgeomfi … … 6 7 SAVE 7 8 9 REAL :: nan ! initialized to NaN with adequate compiler options 10 8 11 REAL, PARAMETER :: pi=2*ASIN(1.), solarcst=1370., stephan=5.67e-08, & 9 12 ps_rad=1.e5, height_scale=10000., ref_temp=285., & 10 capcal_nosoil=1e5 11 12 REAL, ALLOCATABLE :: tsurf(:),tsoil(:,:),rnatur(:), & 13 capcal(:),fluxgrd(:), & 14 dtrad(:,:),fluxrad(:), & 15 q2(:,:),q2l(:,:), & 16 albedo(:),emissiv(:),z0(:),inertie(:) 17 !$OMP THREADPRIVATE( tsurf,tsoil,rnatur) 18 !$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad) 19 !$OMP THREADPRIVATE( q2,q2l) 20 !$OMP THREADPRIVATE( albedo,emissiv,z0,inertie) 13 capcal_nosoil=1e5, tsoil_init=300. 14 15 ! internal state, written to / read from disk at checkpoint / restart 16 REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:), capcal(:), fluxgrd(:) 17 !$OMP THREADPRIVATE(tsurf, tsoil, capcal, fluxgrd) 18 19 ! precomputed arrays 20 REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:) 21 !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie) 21 22 22 23 INTEGER :: icount … … 74 75 75 76 ! Local variables : 76 REAL, DIMENSION(ngrid) :: mu0 ,fract77 REAL, DIMENSION(ngrid) :: mu0 77 78 INTEGER :: j,l,ig,nlevel,igout 78 79 ! … … 88 89 REAL zdum3(ngrid,nlayer) 89 90 REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer) 90 REAL zfluxsw(ngrid),zfluxlw(ngrid) 91 92 print*,'OK DANS PHYPARAM'93 print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid)94 print*,'nlayer',nlayer95 91 REAL zfluxsw(ngrid),zfluxlw(ngrid), fluxrad(ngrid) 92 93 WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid) 94 WRITELOG(*,*) 'nlayer',nlayer 95 LOG_DBG('phyparam') 96 96 97 IF (ngrid.NE.ngridmax) THEN 97 98 PRINT*,'STOP in inifis' … … 113 114 ! zday_last=rjourvrai 114 115 zday_last=zday-ptimestep/unjours 115 CALL alloc_phyparam(ngrid, nlayer, igout) 116 117 ! print*,'OK PHYPARAM 1 ' 118 IF(callsoil) THEN 119 CALL soil(ngrid,nsoilmx,firstcall,inertie, & 120 & ptimestep,tsurf,tsoil,capcal,fluxgrd) 121 ! NB : this call to soil also performs some calculations, see surface.F90 122 ELSE 123 PRINT*,'WARNING!!! Thermal conduction in the soil turned off' 124 DO ig=1,ngrid 125 capcal(ig) = capcal_nosoil 126 fluxgrd(ig) = 0. 127 ENDDO 128 ENDIF 116 CALL alloc(ngrid, nlayer) 117 CALL precompute 118 CALL coldstart(ngrid, ptimestep) 129 119 ENDIF 130 120 … … 172 162 173 163 !----------------------------------------------------------------------- 174 ! 2. C alcul of the radiative tendencies :164 ! 2. Compute radiative tendencies : 175 165 ! --------------------------------------- 176 166 177 167 IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, & 178 gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &179 pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, mu0)168 gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, & 169 pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0) 180 170 181 171 !----------------------------------------------------------------------- … … 184 174 ! 185 175 IF(calldifv) THEN 186 176 187 177 DO ig=1,ngrid 188 178 zflubid(ig)=fluxrad(ig)+fluxgrd(ig) … … 198 188 enddo 199 189 200 CALL vdif(ngrid,nlayer,zday, 201 & ptimestep,capcal,z0, 202 & pplay,pplev,zzlay,zzlev, 203 & pu,pv,zh,tsurf,emissiv, 204 & zdum1,zdum2,zdum3,zflubid, 205 & zdufr,zdvfr,zdhfr,zdtsrfr, q2,q2l,&190 CALL vdif(ngrid,nlayer,zday, & 191 & ptimestep,capcal,z0, & 192 & pplay,pplev,zzlay,zzlev, & 193 & pu,pv,zh,tsurf,emissiv, & 194 & zdum1,zdum2,zdum3,zflubid, & 195 & zdufr,zdvfr,zdhfr,zdtsrfr, & 206 196 & lverbose) 207 197 … … 275 265 & ptimestep,tsurf,tsoil,capcal,fluxgrd) 276 266 IF(lverbose) THEN 277 PRINT*,'Surface Heat capacity,conduction Flux, Ts, dTs, dt' 278 PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout), & 279 & zdtsrf(igout),ptimestep 267 WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt' 268 WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout), & 269 & zdtsrf(igout), ptimestep 270 LOG_DBG('phyparam') 280 271 ENDIF 281 272 ENDIF … … 284 275 ! sorties: 285 276 ! -------- 286 287 print*,'zday, zday_last ',zday,zday_last,icount 277 278 WRITELOG(*,*) 'zday, zday_last ',zday,zday_last,icount 279 LOG_DBG('phyparam') 280 288 281 if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then 289 print*,'zday, zday_last SORTIE ',zday,zday_last 282 WRITELOG(*,*) 'zday, zday_last SORTIE ', zday, zday_last 283 LOG_INFO('phyparam') 284 290 285 zday_last=zday 291 286 ! Ecriture/extension de la coordonnee temps … … 314 309 call writefield('mu0','mu0',' ',mu0) 315 310 call writefield('alb','alb',' ',albedo) 316 call writefield('fract','fract',' ',fract)317 311 call writefield('ps','Surface pressure','Pa',pplev(:,1)) 318 312 call writefield('slp','Sea level pressure','Pa',zpmer) … … 326 320 SUBROUTINE radiative_tendencies(ngrid, igout, nlayer, & 327 321 gmtime, zdtime, zday, pplev, pplay, pt, & 328 pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, mu0)322 pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0) 329 323 USE planet 330 324 USE phys_const, ONLY : planet_rad, unjours … … 339 333 & pt(ngrid, nlayer+1) 340 334 REAL, INTENT(INOUT) :: pdt(ngrid,nlayer) 341 REAL, INTENT(OUT) :: zdtlw(ngrid,nlayer), zfluxlw(ngrid), & 342 zdtsw(ngrid,nlayer), zfluxsw(ngrid), mu0(ngrid) 343 REAL, DIMENSION(ngrid) :: fract 335 REAL, INTENT(OUT) :: zdtlw(ngrid,nlayer), & ! long-wave temperature tendency 336 & zfluxlw(ngrid), & ! long-wave flux at surface 337 & zdtsw(ngrid,nlayer), & ! short-wave temperature tendency 338 & zfluxsw(ngrid), & ! short-wave flux at surface 339 & fluxrad(ngrid), & ! surface flux 340 mu0(ngrid) ! ?? 341 ! local variables 342 REAL :: fract(ngrid), dtrad(ngrid, nlayer) 344 343 REAL :: zls, zinsol, tsurf2, ztim1,ztim2,ztim3, dist_sol, declin 345 344 REAL :: zplanck(ngrid) … … 362 361 ELSE 363 362 CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract) 364 print*,'ZENANG '365 363 ENDIF 366 364 367 365 IF(lverbose) THEN 368 PRINT*,'day, declin, sinlon,coslon,sinlat,coslat' 369 PRINT*,zday, declin, & 370 & sinlon(igout),coslon(igout), & 371 & sinlat(igout),coslat(igout) 366 WRITELOG(*,*) 'day, declin, sinlon,coslon,sinlat,coslat' 367 WRITELOG(*,*) zday, declin, & 368 & sinlon(igout),coslon(igout), & 369 & sinlat(igout),coslat(igout) 370 LOG_DBG('radiative_tendencies') 372 371 ENDIF 373 372 ELSE 374 print*,'declin,ngrid,planet_rad',declin,ngrid,planet_rad 373 WRITELOG(*,*) 'declin,ngrid,planet_rad', declin, ngrid, planet_rad 374 LOG_DBG('radiative_tendencies') 375 375 CALL mucorr(ngrid,declin,lati,mu0,fract,height_scale,planet_rad) 376 376 ENDIF … … 414 414 415 415 IF(lverbose) THEN 416 PRINT*,'Diagnostics for radiation'417 PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'418 PRINT*,albedo(igout),emissiv(igout),mu0(igout), &419 & 420 & 421 PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'422 PRINT*,'unjours',unjours416 WRITELOG(*,*) 'Diagnostics for radiation' 417 WRITELOG(*,*) 'albedo, emissiv, mu0,fract,Frad,Planck' 418 WRITELOG(*,*) albedo(igout),emissiv(igout),mu0(igout), & 419 & fract(igout), & 420 & fluxrad(igout),zplanck(igout) 421 WRITELOG(*,*) 'Tlay Play Plev dT/dt SW dT/dt LW (K/day)' 422 WRITELOG(*,*) 'unjours',unjours 423 423 DO l=1,nlayer 424 WRITE (*,'(3f15.5,2e15.2)') pt(igout,l),&424 WRITELOG(*,'(3f15.5,2e15.2)') pt(igout,l), & 425 425 & pplay(igout,l),pplev(igout,l), & 426 426 & zdtsw(igout,l),zdtlw(igout,l) 427 427 ENDDO 428 LOG_DBG('radiative_tendencies') 428 429 ENDIF 429 430 430 431 END SUBROUTINE radiative_tendencies 431 432 432 SUBROUTINE alloc_phyparam(ngrid, nlayer, igout) 433 SUBROUTINE alloc(ngrid, nlayer) 434 USE astronomy, ONLY : iniorbit 435 USE surface, ONLY : zc,zd 436 INTEGER, INTENT(IN) :: ngrid, nlayer 437 LOGICAL, PARAMETER :: firstcall=.TRUE. 438 ! allocate arrays for internal state 439 ALLOCATE(tsurf(ngrid)) 440 ALLOCATE(tsoil(ngrid,nsoilmx)) 441 ! we could avoid the arrays below with a different implementation of surface / radiation / turbulence coupling 442 ALLOCATE(capcal(ngrid),fluxgrd(ngrid)) 443 ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx)) 444 ! allocate precomputed arrays 445 ALLOCATE(rnatur(ngrid)) 446 ALLOCATE(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid)) 447 CALL iniorbit 448 END SUBROUTINE alloc 449 450 SUBROUTINE precompute 433 451 USE surface 434 USE astronomy, ONLY : iniorbit 435 INTEGER, INTENT(IN) :: ngrid, nlayer, igout 436 LOGICAL, PARAMETER :: firstcall=.TRUE. 437 438 print*,'AKk',ngrid,nsoilmx 439 allocate(tsurf(ngrid)) 440 print*,'AKa' 441 allocate (tsoil(ngrid,nsoilmx)) 442 print*,'AKb' 443 allocate (rnatur(ngrid)) 444 print*,'AK2' 445 allocate(capcal(ngrid),fluxgrd(ngrid)) 446 print*,'AK3' 447 allocate(dtrad(ngrid,nlayer),fluxrad(ngrid)) 448 print*,'AK4' 449 allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1)) 450 print*,'AK5' 451 allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid)) 452 print*,'AK6' 453 454 455 PRINT*,'FIRSTCALL ' 456 457 rnatur=1. 458 emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter 459 inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter 460 albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter 461 z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter 462 q2=1.e-10 463 q2l=1.e-10 464 tsurf=300. 465 tsoil=300. 466 print*,tsoil(igout,nsoilmx/2+2) 467 print*,'TS ',tsurf(igout),tsoil(igout,5) 468 CALL iniorbit 469 470 if (.not.callrad) fluxrad(:)=0. 471 452 ! precompute time-independent arrays 453 rnatur(:) = 1. 454 emissiv(:) = (1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter 455 inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter 456 albedo(:) = (1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter 457 z0(:) = (1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter 458 END SUBROUTINE precompute 459 460 SUBROUTINE coldstart(ngrid, ptimestep) 461 ! create internal state to start a run without a restart file 462 USE surface, ONLY : soil 463 INTEGER, INTENT(IN) :: ngrid 464 REAL, INTENT(iN) :: ptimestep 465 tsurf(:) = tsoil_init 466 tsoil(:,:) = tsoil_init 472 467 icount=0 473 474 END SUBROUTINE alloc_phyparam 468 IF(callsoil) THEN 469 ! initializes zc, zd, capcal, fluxgrd 470 CALL soil(ngrid,nsoilmx,.TRUE.,inertie, & 471 & ptimestep,tsurf,tsoil,capcal,fluxgrd) 472 ELSE 473 WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off' 474 LOG_WARN('coldstart') 475 capcal(:) = capcal_nosoil 476 fluxgrd(:) = 0. 477 ENDIF 478 END SUBROUTINE coldstart 475 479 476 480 END MODULE phyparam_mod -
dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90
r4210 r4228 53 53 ! ---------- 54 54 ! 55 INTEGER ngrid,nlayer 56 REAL albedo(ngrid),coefvis 57 REAL pmu(ngrid),pfract(ngrid) 58 REAL plevel(ngrid,nlayer+1),ps_rad 59 REAL psolarf0 60 REAL fsrfvis(ngrid),dtsw(ngrid,nlayer) 61 LOGICAL lwrite,ldiurn 62 ! 63 ! variables locales: 64 ! ------------------ 65 ! 55 INTEGER, INTENT(IN) :: ngrid,nlayer 56 LOGICAL, INTENT(IN) :: lwrite,ldiurn 57 REAL, INTENT(IN) :: albedo(ngrid), coefvis 58 REAL, INTENT(IN) :: pmu(ngrid), pfract(ngrid) 59 REAL, INTENT(IN) :: plevel(ngrid,nlayer+1), ps_rad 60 REAL, INTENT(IN) :: psolarf0 61 REAL, INTENT(OUT) :: fsrfvis(ngrid),dtsw(ngrid,nlayer) 66 62 67 63 REAL zalb(ngrid),zmu(ngrid),zfract(ngrid) … … 80 76 ! ------------------- 81 77 82 83 78 nlevel=nlayer+1 84 79 85 80 !----------------------------------------------------------------------- 86 81 ! Definitions des tableaux locaux pour les points ensoleilles: -
dynamico_lmdz/simple_physics/phyparam/physics/surface.F90
r4210 r4228 13 13 & alb_mer,alb_ter,emi_mer,emi_ter 14 14 15 ! local saved variables: 16 ! ---------------------- 15 ! precomputed variables 17 16 REAL :: lambda 18 REAL,ALLOCATABLE :: dz1(:),dz2(:),zc(:,:),zd(:,:) 19 !$OMP THREADPRIVATE(dz1,dz2,zc,zd,lambda) 20 21 PUBLIC :: soil 17 REAL,ALLOCATABLE :: dz1(:),dz2(:) 18 !$OMP THREADPRIVATE(dz1,dz2,zc,lambda) 19 20 ! internal state variables needed for checkpoint/restart 21 REAL, ALLOCATABLE :: zc(:,:),zd(:,:) 22 !$OMP THREADPRIVATE(zc,zd) 23 24 PUBLIC :: soil, zc, zd 22 25 23 26 CONTAINS … … 40 43 41 44 ALLOCATE(dz1(nsoil),dz2(nsoil)) 42 ALLOCATE(zc(ngrid,nsoil),zd(ngrid,nsoil))43 45 44 46 min_period=20000. … … 59 61 ENDDO 60 62 lambda=fz(.5)*dz1(1) 63 61 64 WRITELOG(*,*) 'full layers, intermediate layers (secoonds)' 62 65 DO jk=1,nsoil … … 67 70 & fz(rk)*fz(rk)*pi 68 71 ENDDO 69 70 72 LOG_INFO('init_soil') 73 71 74 END SUBROUTINE init_soil 72 75 … … 134 137 135 138 INTEGER ig,jk 136 REAL za(ngrid),zb(ngrid)137 139 REAL zdz2(nsoil),z1(ngrid) 138 140 … … 154 156 DO jk=1,nsoil-1 155 157 DO ig=1,ngrid 156 ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) 158 ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) 157 159 ENDDO 158 160 ENDDO -
dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90
r4225 r4228 101 101 ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pkv,pkh) 102 102 ! FIXME : pkh := pkv 103 USE planet104 103 INTEGER, INTENT(IN) :: ngrid,nlay 105 104 REAL, INTENT(IN) :: ptimestep, pg, & … … 160 159 pu,pv,ph,ptsrf,pemis, & 161 160 pdufi,pdvfi,pdhfi,pfluxsrf, & 162 pdudif,pdvdif,pdhdif,pdtsrf, pq2,pq2l,&161 pdudif,pdvdif,pdhdif,pdtsrf, & 163 162 lwrite) 164 163 USE phys_const … … 195 194 REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay) 196 195 REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid) 197 REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)198 196 LOGICAL lwrite 199 197 !
Note: See TracChangeset
for help on using the changeset viewer.