Changeset 4212 for dynamico_lmdz
- Timestamp:
- Jan 7, 2020, 12:46:42 PM (5 years ago)
- Location:
- dynamico_lmdz/simple_physics/phyparam
- Files:
-
- 1 deleted
- 3 edited
- 3 moved
-
dynphy_lonlat/phyaqua.F (moved) (moved from dynamico_lmdz/simple_physics/phyparam/param/phyaqua.F)
-
param/callkeys.F90 (moved) (moved from dynamico_lmdz/simple_physics/phyparam/param/callkeys.h) (1 diff)
-
param/clesphys.h (deleted)
-
param/comgeomfi.F90 (modified) (2 diffs)
-
param/iniphyparam.F (modified) (2 diffs)
-
param/phyparam.F90 (moved) (moved from dynamico_lmdz/simple_physics/phyparam/param/phyparam.F) (5 diffs)
-
param/physiq_mod.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/phyparam/param/callkeys.F90
r4196 r4212 1 COMMON/callkeys/callrad,calldifv,calladj,callcond,callsoil, 2 s season,diurnal,lverbose,iradia,period_sort 3 LOGICAL callrad,calldifv,calladj,callcond,callsoil, 4 s season,diurnal,lverbose 5 INTEGER iradia 6 real period_sort 1 MODULE callkeys 2 IMPLICIT NONE 3 SAVE 4 5 LOGICAL callrad,calldifv,calladj,callcond,callsoil,season,diurnal,lverbose 6 INTEGER iradia 7 REAL period_sort 8 9 END MODULE callkeys -
dynamico_lmdz/simple_physics/phyparam/param/comgeomfi.F90
r4176 r4212 1 module comgeomfi 1 MODULE comgeomfi 2 IMPLICIT NONE 3 SAVE 2 4 3 real,save,allocatable :: long(:) 4 real,save,allocatable :: lati(:) 5 real,save,allocatable :: area(:) 6 real,save,allocatable :: sinlon(:) 7 real,save,allocatable :: coslon(:) 8 real,save,allocatable :: sinlat(:) 9 real,save,allocatable :: coslat(:) 10 real,save :: totarea 11 integer,save :: ngridmax,nlayermx,nsoilmx 5 REAL, allocatable :: long(:), lati(:), area(:), & 6 sinlon(:), coslon(:), sinlat(:), coslat(:) 7 REAL :: totarea 8 INTEGER :: ngridmax,nlayermx,nsoilmx 12 9 !$OMP THREADPRIVATE(long,lati,area,sinlon,coslon,sinlat,coslat,totarea) 13 10 !$OMP THREADPRIVATE(ngridmax,nlayermx,nsoilmx) 14 contains15 16 subroutine InitComgeomfi17 USE mod_phys_lmdz_para18 USE dimphy, ONLY : klon,klev19 USE geometry_mod, ONLY : latitude_deg,longitude_deg20 implicit none21 11 12 CONTAINS 13 14 SUBROUTINE InitComgeomfi 15 USE mod_phys_lmdz_para 16 USE dimphy, ONLY : klon,klev 17 USE geometry_mod, ONLY : latitude_deg,longitude_deg 22 18 23 19 print*,'Dans initcomgeomfi ' … … 37 33 allocate(coslat(klon_omp)) 38 34 39 end subroutineInitComgeomfi35 END SUBROUTINE InitComgeomfi 40 36 41 end modulecomgeomfi37 END MODULE comgeomfi -
dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F
r4208 r4212 8 8 USE mod_grid_phy_lmdz 9 9 USE mod_phys_lmdz_para 10 USE callkeys 10 11 use comgeomfi 11 12 use comsaison … … 52 53 53 54 54 #include "callkeys.h"55 55 #include "iniprint.h" 56 56 -
dynamico_lmdz/simple_physics/phyparam/param/phyparam.F90
r4208 r4212 1 SUBROUTINE phyparam(ngrid,nlayer,nq, 2 s firstcall,lastcall, 3 s rjourvrai,gmtime,ptimestep, 4 s pplev,pplay,pphi,pphis,presnivs, 5 s pu,pv,pt,pq, 6 s pw, 7 s pdu,pdv,pdt,pdq,pdpsrf) 8 9 USE comsaison 10 USE dimphy 11 USE comgeomfi 12 USE phys_const 13 USE planet 14 USE astronomy 15 USE turbulence, ONLY : vdif 16 USE convection, ONLY : convadj 17 USE solar, ONLY : solang, zenang, mucorr 18 USE radiative_sw, ONLY : sw 19 USE radiative_lw, ONLY : lw 20 USE surface 21 c 22 IMPLICIT NONE 23 c======================================================================= 24 c 25 c subject: 26 c -------- 27 c 28 c Organisation of the physical parametrisations of the LMD 29 c 20 parameters GCM for planetary atmospheres. 30 c It includes: 31 c raditive transfer (long and shortwave) for CO2 and dust. 32 c vertical turbulent mixing 33 c convective adjsutment 34 c 35 c author: Frederic Hourdin 15 / 10 /93 36 c ------- 37 c 38 c arguments: 39 c ---------- 40 c 41 c input: 42 c ------ 43 c 44 c ngrid Size of the horizontal grid. 45 c All internal loops are performed on that grid. 46 c nlayer Number of vertical layers. 47 c nq Number of advected fields 48 c firstcall True at the first call 49 c lastcall True at the last call 50 c rjourvrai Number of days counted from the North. Spring 51 c equinoxe. 52 c gmtime hour (s) 53 c ptimestep timestep (s) 54 c pplay(ngrid,nlayer+1) Pressure at the middle of the layers (Pa) 55 c pplev(ngrid,nlayer+1) intermediate pressure levels (pa) 56 c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2) 57 c pu(ngrid,nlayer) u component of the wind (ms-1) 58 c pv(ngrid,nlayer) v component of the wind (ms-1) 59 c pt(ngrid,nlayer) Temperature (K) 60 c pq(ngrid,nlayer,nq) Advected fields 61 c pudyn(ngrid,nlayer) \ 62 c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the 63 c ptdyn(ngrid,nlayer) / corresponding variables 64 c pqdyn(ngrid,nlayer,nq) / 65 c pw(ngrid,?) vertical velocity 66 c 67 c output: 68 c ------- 69 c 70 c pdu(ngrid,nlayermx) \ 71 c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding 72 c pdt(ngrid,nlayermx) / variables due to physical processes. 73 c pdq(ngrid,nlayermx) / 74 c pdpsrf(ngrid) / 75 c 76 c======================================================================= 77 c 78 c----------------------------------------------------------------------- 79 c 80 c 0. Declarations : 81 c ------------------ 1 MODULE phyparam_mod 2 IMPLICIT NONE 3 4 CONTAINS 5 6 SUBROUTINE phyparam(ngrid,nlayer,nq, & 7 & firstcall,lastcall, & 8 & rjourvrai,gmtime,ptimestep, & 9 & pplev,pplay,pphi,pphis,presnivs, & 10 & pu,pv,pt,pq, & 11 & pw, & 12 & pdu,pdv,pdt,pdq,pdpsrf) 13 14 USE callkeys 15 USE comsaison 16 USE dimphy 17 USE comgeomfi 18 USE phys_const 19 USE planet 20 USE astronomy 21 USE turbulence, ONLY : vdif 22 USE convection, ONLY : convadj 23 USE solar, ONLY : solang, zenang, mucorr 24 USE radiative_sw, ONLY : sw 25 USE radiative_lw, ONLY : lw 26 USE surface 27 ! 28 !======================================================================= 29 ! 30 ! subject: 31 ! -------- 32 ! 33 ! Organisation of the physical parametrisations of the LMD 34 ! 20 parameters GCM for planetary atmospheres. 35 ! It includes: 36 ! raditive transfer (long and shortwave) for CO2 and dust. 37 ! vertical turbulent mixing 38 ! convective adjsutment 39 ! 40 ! author: Frederic Hourdin 15 / 10 /93 41 ! ------- 42 ! 43 ! arguments: 44 ! ---------- 45 ! 46 ! input: 47 ! ------ 48 ! 49 ! ngrid Size of the horizontal grid. 50 ! All internal loops are performed on that grid. 51 ! nlayer Number of vertical layers. 52 ! nq Number of advected fields 53 ! firstcall True at the first call 54 ! lastcall True at the last call 55 ! rjourvrai Number of days counted from the North. Spring 56 ! equinoxe. 57 ! gmtime hour (s) 58 ! ptimestep timestep (s) 59 ! pplay(ngrid,nlayer+1) Pressure at the middle of the layers (Pa) 60 ! pplev(ngrid,nlayer+1) intermediate pressure levels (pa) 61 ! pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2) 62 ! pu(ngrid,nlayer) u component of the wind (ms-1) 63 ! pv(ngrid,nlayer) v component of the wind (ms-1) 64 ! pt(ngrid,nlayer) Temperature (K) 65 ! pq(ngrid,nlayer,nq) Advected fields 66 ! pudyn(ngrid,nlayer) \ 67 ! pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the 68 ! ptdyn(ngrid,nlayer) / corresponding variables 69 ! pqdyn(ngrid,nlayer,nq) / 70 ! pw(ngrid,?) vertical velocity 71 ! 72 ! output: 73 ! ------- 74 ! 75 ! pdu(ngrid,nlayermx) \ 76 ! pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding 77 ! pdt(ngrid,nlayermx) / variables due to physical processes. 78 ! pdq(ngrid,nlayermx) / 79 ! pdpsrf(ngrid) / 80 ! 81 !======================================================================= 82 ! 83 !----------------------------------------------------------------------- 84 ! 85 ! 0. Declarations : 86 ! ------------------ 82 87 83 88 #include "dimensions.h" 84 #include "callkeys.h" 85 c#include "surface.h" 86 87 c Arguments : 88 c ----------- 89 90 c inputs: 91 c ------- 92 INTEGER ngrid,nlayer,nq 93 94 REAL ptimestep 95 real zdtime 96 REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer) 97 REAL pphi(ngrid,nlayer) 98 REAL pphis(ngrid) 99 REAL pu(ngrid,nlayer),pv(ngrid,nlayer) 100 REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq) 101 REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer) 102 103 c dynamial tendencies 104 REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq) 105 REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer) 106 REAL pw(ngrid,nlayer) 107 108 c Time 109 real rjourvrai 110 REAL gmtime 111 112 c outputs: 113 c -------- 114 115 c physical tendencies 116 REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq) 117 REAL pdpsrf(ngrid) 118 LOGICAL firstcall,lastcall 119 120 c Local variables : 121 c ----------------- 122 123 INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,isoil,iq 124 INTEGER*4 day_ini 125 c 126 REAL,DIMENSION(ngrid) :: mu0,fract 127 REAL zday 128 REAL zh(ngrid,nlayer),z1,z2 129 REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer) 130 REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer) 131 REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid) 132 REAL zflubid(ngrid),zpmer(ngrid) 133 REAL zplanck(ngrid),zpopsk(ngrid,nlayer) 134 REAL zdum1(ngrid,nlayer) 135 REAL zdum2(ngrid,nlayer) 136 REAL zdum3(ngrid,nlayer) 137 REAL ztim1,ztim2,ztim3 138 REAL zls,zinsol 139 REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer) 140 REAL zfluxsw(ngrid),zfluxlw(ngrid) 141 character*2 str2 142 REAL factq(nq),tauq(nq) 143 character*3 nomq 144 145 c Local saved variables: 146 c ---------------------- 147 148 INTEGER, SAVE :: icount 149 real, SAVE :: zday_last 150 !$OMP THREADPRIVATE( icount,zday_last) 151 152 REAL zps_av 153 154 real,allocatable,SAVE :: tsurf(:),tsoil(:,:),rnatur(:) 155 real,allocatable,SAVE :: capcal(:),fluxgrd(:) 156 real,allocatable,SAVE :: dtrad(:,:),fluxrad(:) 157 real,allocatable,SAVE :: q2(:,:),q2l(:,:) 158 real,allocatable,SAVE :: albedo(:),emissiv(:),z0(:),inertie(:) 159 real,SAVE :: solarcst=1370. 160 real,SAVE :: stephan=5.67e-08 161 162 !$OMP THREADPRIVATE(tsurf,tsoil,rnatur) 163 !$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad) 164 !$OMP THREADPRIVATE( q2,q2l) 165 !$OMP THREADPRIVATE( albedo,emissiv,solarcst,z0,inertie) 166 !$OMP THREADPRIVATE( stephan) 167 168 169 EXTERNAL ismin,ismax 170 171 172 INTEGER longcles 173 PARAMETER ( longcles = 20 ) 174 REAL clesphy0( longcles ) 175 REAL presnivs(nlayer) 176 177 178 179 print*,'OK DANS PHYPARAM' 180 181 c----------------------------------------------------------------------- 182 c 1. Initialisations : 183 c -------------------- 184 185 ! call initial0(ngrid*nlayermx*nqmx,pdq) 186 nlevel=nlayer+1 187 188 ! print*,'OK ',nlevel 189 190 igout=ngrid/2+1 191 ! print*,'OK PHYPARAM ',ngrid,igout 192 193 194 zday=rjourvrai+gmtime 195 196 ! print*,'OK PHYPARAM 0A nq ',nq 197 tauq(1)=1800. 198 tauq(2)=10800. 199 tauq(3)=86400. 200 tauq(4)=864000. 201 ! print*,'OK PHYPARAM 0 B nq ',nq 202 factq(1:4)=(1.-exp(-ptimestep/tauq(1:4)))/ptimestep 203 204 ! print*,'OK PHYPARAM 0 ' 205 print*,'nq ',nq 206 print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid) 207 print*,'nlayer',nlayer 208 print*,'size pdq ',ngrid*nlayer*4,ngrid*nlayer*nq, 209 s size(pdq),size(lati),size(pq),size(factq) 210 do iq=1,4 211 do l=1,nlayer 212 pdq(1:ngrid,l,iq)= 213 & (1.+sin(lati(1:ngrid))-pq(1:ngrid,l,iq))*factq(iq) 214 enddo 215 enddo 216 217 IF(firstcall) THEN 218 219 print*,'AKk',ngrid,nsoilmx 220 allocate(tsurf(ngrid)) 221 print*,'AKa' 222 allocate (tsoil(ngrid,nsoilmx)) 223 print*,'AKb' 224 allocate (rnatur(ngrid)) 225 print*,'AK2' 226 allocate(capcal(ngrid),fluxgrd(ngrid)) 227 print*,'AK3' 228 allocate(dtrad(ngrid,nlayer),fluxrad(ngrid)) 229 print*,'AK4' 230 allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1)) 231 print*,'AK5' 232 allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid)) 233 print*,'AK6' 234 235 236 do l=1,nlayer 237 pdq(:,l,5)=1.+sin(lati(:))/ptimestep 238 enddo 239 PRINT*,'FIRSTCALL ' 240 241 ! zday_last=rjourvrai 242 zday_last=zday-ptimestep/unjours 243 c CALL readfi(ngrid,nlayer,nsoilmx,ldrs, 244 c . day_ini,gmtime,albedo,inertie,emissiv,z0,rnatur, 245 c . q2,q2l,tsurf,tsoil) 246 rnatur=1. 247 emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter 248 inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter 249 albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter 250 z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter 251 q2=1.e-10 252 q2l=1.e-10 253 tsurf=300. 254 tsoil=300. 255 print*,tsoil(ngrid/2+1,nsoilmx/2+2) 256 print*,'TS ',tsurf(igout),tsoil(igout,5) 257 CALL iniorbit 258 259 if (.not.callrad) then 260 do ig=1,ngrid 261 fluxrad(ig)=0. 262 enddo 263 endif 264 265 ! print*,'OK PHYPARAM 1 ' 266 IF(callsoil) THEN 267 CALL soil(ngrid,nsoilmx,firstcall,inertie, 268 s ptimestep,tsurf,tsoil,capcal,fluxgrd) 269 ELSE 270 PRINT*,'WARNING!!! Thermal conduction in the soil 271 s turned off' 272 DO ig=1,ngrid 273 capcal(ig)=1.e5 274 fluxgrd(ig)=0. 275 ENDDO 276 ENDIF 277 c CALL inifrict(ptimestep) 278 icount=0 279 280 PRINT*,'FIRSTCALL B ' 281 print*,'INIIO avant iophys_ini ' 282 call iophys_ini('phys.nc ',presnivs) 283 284 ENDIF 285 icount=icount+1 286 287 PRINT*,'FIRSTCALL AP ' 288 IF (ngrid.NE.ngridmax) THEN 289 PRINT*,'STOP in inifis' 290 PRINT*,'Probleme de dimenesions :' 291 PRINT*,'ngrid = ',ngrid 292 PRINT*,'ngridmax = ',ngridmax 293 STOP 294 ENDIF 295 296 DO l=1,nlayer 297 DO ig=1,ngrid 298 pdv(ig,l)=0. 299 pdu(ig,l)=0. 300 pdt(ig,l)=0. 301 ENDDO 302 ENDDO 303 304 DO ig=1,ngrid 305 pdpsrf(ig)=0. 306 zflubid(ig)=0. 307 zdtsrf(ig)=0. 308 ENDDO 309 310 c----------------------------------------------------------------------- 311 c calcul du geopotentiel aux niveaux intercouches 312 c ponderation des altitudes au niveau des couches en dp/p 313 314 DO l=1,nlayer 315 DO ig=1,ngrid 316 zzlay(ig,l)=pphi(ig,l)/g 317 ENDDO 318 ENDDO 319 DO ig=1,ngrid 320 zzlev(ig,1)=0. 321 ENDDO 322 DO l=2,nlayer 323 DO ig=1,ngrid 324 z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) 325 z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) 326 zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) 327 ENDDO 328 ENDDO 329 330 c----------------------------------------------------------------------- 331 c Transformation de la temperature en temperature potentielle 332 DO l=1,nlayer 333 DO ig=1,ngrid 334 zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp 335 ENDDO 336 ENDDO 337 DO l=1,nlayer 338 DO ig=1,ngrid 339 zh(ig,l)=pt(ig,l)/zpopsk(ig,l) 340 ENDDO 341 ENDDO 342 343 c----------------------------------------------------------------------- 344 c 2. Calcul of the radiative tendencies : 345 c --------------------------------------- 346 347 ! print*,'callrad0' 348 IF(callrad) THEN 349 ! print*,'callrad' 350 351 c WARNING !!! on calcule le ray a chaque appel 352 c IF( MOD(icount,iradia).EQ.0) THEN 353 354 CALL solarlong(zday,zls) 355 CALL orbite(zls,dist_sol,declin) 356 c declin=0. 357 ! print*,'ATTENTIOn : pdeclin = 0',' L_s=',zls 358 print*,'diurnal=',diurnal 359 IF(diurnal) THEN 360 if ( 1.eq.1 ) then 361 ztim1=SIN(declin) 362 ztim2=COS(declin)*COS(2.*pi*(zday-.5)) 363 ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) 364 c call dump2d(iim,jjm-1,lati(2),'LATI 0 ') 365 c call dump2d(iim,jjm-1,long(2),'LONG 0 ') 366 c call dump2d(iim,jjm-1,sinlon(2),'sinlon0 ') 367 c call dump2d(iim,jjm-1,coslon(2),'coslon0 ') 368 c call dump2d(iim,jjm-1,sinlat(2),'sinlat ') 369 c call dump2d(iim,jjm-1,coslat(2),'coslat ') 370 CALL solang(ngrid,sinlon,coslon,sinlat,coslat, 371 s ztim1,ztim2,ztim3, 372 s mu0,fract) 373 else 374 zdtime=ptimestep*float(iradia) 375 CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract) 376 print*,'ZENANG ' 377 endif 378 379 IF(lverbose) THEN 380 PRINT*,'day, declin, sinlon,coslon,sinlat,coslat' 381 PRINT*,zday, declin, 382 s sinlon(igout),coslon(igout), 383 s sinlat(igout),coslat(igout) 384 ENDIF 385 ELSE 386 print*,'declin,ngrid,rad',declin,ngrid,rad 387 388 c call dump2d(iim,jjm-1,lati(2),'LATI ') 389 CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad) 390 ENDIF 391 c call dump2d(iim,jjm-1,fract(2),'FRACT A ') 392 c call dump2d(iim,jjm-1,mu0(2),'MU0 A ') 393 394 395 c 2.2 Calcul of the radiative tendencies and fluxes: 396 c -------------------------------------------------- 397 398 c 2.1.2 levels 399 400 zinsol=solarcst/(dist_sol*dist_sol) 401 print*,iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer 402 print*,'iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer' 403 c call dump2d(iim,jjm-1,albedo(2),'ALBEDO ') 404 c call dump2d(iim,jjm-1,mu0(2),'MU0 ') 405 c call dump2d(iim,jjm-1,fract(2),'FRACT ') 406 c call dump2d(iim,jjm-1,lati(2),'LATI ') 407 zps_av=1.e5 408 if (firstcall) then 409 print*,'WARNING ps_rad impose' 410 endif 411 CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, 412 $ pplev,zps_av, 413 $ mu0,fract,zinsol, 414 $ zfluxsw,zdtsw, 415 $ lverbose) 416 c call dump2d(iim,jjm-1,zfluxsw(2),'SWS 1 ') 417 c stop 418 419 CALL lw(ngrid,nlayer,coefir,emissiv, 420 $ pplev,zps_av,tsurf,pt, 421 $ zfluxlw,zdtlw, 422 $ lverbose) 423 424 425 c 2.4 total flux and tendencies: 426 c ------------------------------ 427 428 c 2.4.1 fluxes 429 430 DO ig=1,ngrid 431 fluxrad(ig)=emissiv(ig)*zfluxlw(ig) 432 $ +zfluxsw(ig)*(1.-albedo(ig)) 433 zplanck(ig)=tsurf(ig)*tsurf(ig) 434 zplanck(ig)=emissiv(ig)* 435 $ stephan*zplanck(ig)*zplanck(ig) 436 fluxrad(ig)=fluxrad(ig)-zplanck(ig) 437 ENDDO 438 439 c 2.4.2 temperature tendencies 440 441 DO l=1,nlayer 442 DO ig=1,ngrid 443 dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) 444 ENDDO 445 ENDDO 446 447 c ENDIF 448 449 c 2.5 Transformation of the radiative tendencies: 450 c ----------------------------------------------- 451 452 DO l=1,nlayer 453 DO ig=1,ngrid 454 pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) 455 ENDDO 456 ENDDO 457 458 IF(lverbose) THEN 459 PRINT*,'Diagnotique for the radaition' 460 PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck' 461 PRINT*,albedo(igout),emissiv(igout),mu0(igout), 462 s fract(igout), 463 s fluxrad(igout),zplanck(igout) 464 PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)' 465 PRINT*,'unjours',unjours 466 DO l=1,nlayer 467 WRITE(*,'(3f15.5,2e15.2)') pt(igout,l), 468 s pplay(igout,l),pplev(igout,l), 469 s zdtsw(igout,l),zdtlw(igout,l) 470 ENDDO 471 ENDIF 472 473 474 ENDIF 475 476 c----------------------------------------------------------------------- 477 c 3. Vertical diffusion (turbulent mixing): 478 c ----------------------------------------- 479 c 480 IF(calldifv) THEN 481 482 DO ig=1,ngrid 483 zflubid(ig)=fluxrad(ig)+fluxgrd(ig) 484 ENDDO 485 486 zdum1(:,:)=0. 487 zdum2(:,:)=0. 488 489 do l=1,nlayer 490 do ig=1,ngrid 491 zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l) 492 enddo 493 enddo 494 495 CALL vdif(ngrid,nlayer,zday, 496 $ ptimestep,capcal,z0, 497 $ pplay,pplev,zzlay,zzlev, 498 $ pu,pv,zh,tsurf,emissiv, 499 $ zdum1,zdum2,zdum3,zflubid, 500 $ zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l, 501 $ lverbose) 502 c igout=ngrid/2+1 503 c PRINT*,'zdufr zdvfr zdhfr' 504 c DO l=1,nlayer 505 c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l) 506 c ENDDO 507 c CALL difv (ngrid,nlayer, 508 c $ area,lati,capcal, 509 c $ pplev,pphi, 510 c $ pu,pv,zh,tsurf,emissiv, 511 c $ zdum1,zdum2,zdum3,zflubid, 512 c $ zdufr,zdvfr,zdhfr,zdtsrf, 513 c $ .true.) 514 c PRINT*,'zdufr zdvfr zdhfr' 515 c DO l=1,nlayer 516 c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l) 517 c ENDDO 518 c STOP 519 520 DO l=1,nlayer 521 DO ig=1,ngrid 522 pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l) 523 pdu(ig,l)=pdu(ig,l)+zdufr(ig,l) 524 pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l) 525 ENDDO 526 ENDDO 527 528 DO ig=1,ngrid 529 zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig) 530 ENDDO 531 532 ELSE 533 DO ig=1,ngrid 534 zdtsrf(ig)=zdtsrf(ig)+ 535 s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) 536 ENDDO 537 ENDIF 538 c 539 c----------------------------------------------------------------------- 540 c 4. Dry convective adjustment: 541 c ----------------------------- 542 543 IF(calladj) THEN 544 545 DO l=1,nlayer 546 DO ig=1,ngrid 547 zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l) 548 ENDDO 549 ENDDO 550 551 zdufr(:,:)=0. 552 zdvfr(:,:)=0. 553 zdhfr(:,:)=0. 554 555 CALL convadj(ngrid,nlayer,ptimestep, 556 $ pplay,pplev,zpopsk, 557 $ pu,pv,zh, 558 $ pdu,pdv,zdum1, 559 $ zdufr,zdvfr,zdhfr) 560 561 DO l=1,nlayer 562 DO ig=1,ngrid 563 pdu(ig,l)=pdu(ig,l)+zdufr(ig,l) 564 pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l) 565 pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l) 566 ENDDO 567 ENDDO 568 569 ENDIF 570 571 c----------------------------------------------------------------------- 572 c On incremente les tendances physiques de la temperature du sol: 573 c --------------------------------------------------------------- 574 575 DO ig=1,ngrid 576 tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) 577 ENDDO 578 579 WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1) 580 581 c----------------------------------------------------------------------- 582 c soil temperatures: 583 c -------------------- 584 585 IF (callsoil) THEN 586 CALL soil(ngrid,nsoilmx,.false.,inertie, 587 s ptimestep,tsurf,tsoil,capcal,fluxgrd) 588 IF(lverbose) THEN 589 PRINT*,'Surface Heat capacity,conduction Flux, Ts, 590 s dTs, dt' 591 PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout), 592 s zdtsrf(igout),ptimestep 593 ENDIF 594 ENDIF 595 596 c----------------------------------------------------------------------- 597 c sorties: 598 c -------- 599 600 c call dump2d(iim,jjm-1,zfluxsw(2),'SWS 2 ') 601 print*,'zday, zday_last ',zday,zday_last,icount 602 if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then 603 print*,'zday, zday_last SORTIE ',zday,zday_last 89 90 ! Arguments : 91 ! ----------- 92 93 ! inputs: 94 ! ------- 95 INTEGER ngrid,nlayer,nq 96 97 REAL ptimestep 98 real zdtime 99 REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer) 100 REAL pphi(ngrid,nlayer) 101 REAL pphis(ngrid) 102 REAL pu(ngrid,nlayer),pv(ngrid,nlayer) 103 REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq) 104 REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer) 105 106 ! dynamial tendencies 107 REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq) 108 REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer) 109 REAL pw(ngrid,nlayer) 110 111 ! Time 112 real rjourvrai 113 REAL gmtime 114 115 ! outputs: 116 ! -------- 117 118 ! physical tendencies 119 REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq) 120 REAL pdpsrf(ngrid) 121 LOGICAL firstcall,lastcall 122 123 ! Local variables : 124 ! ----------------- 125 126 INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,isoil,iq 127 INTEGER*4 day_ini 128 ! 129 REAL,DIMENSION(ngrid) :: mu0,fract 130 REAL zday 131 REAL zh(ngrid,nlayer),z1,z2 132 REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer) 133 REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer) 134 REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid) 135 REAL zflubid(ngrid),zpmer(ngrid) 136 REAL zplanck(ngrid),zpopsk(ngrid,nlayer) 137 REAL zdum1(ngrid,nlayer) 138 REAL zdum2(ngrid,nlayer) 139 REAL zdum3(ngrid,nlayer) 140 REAL ztim1,ztim2,ztim3 141 REAL zls,zinsol 142 REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer) 143 REAL zfluxsw(ngrid),zfluxlw(ngrid) 144 character*2 str2 145 REAL factq(nq),tauq(nq) 146 character*3 nomq 147 148 ! Local saved variables: 149 ! ---------------------- 150 151 INTEGER, SAVE :: icount 152 real, SAVE :: zday_last 153 !$OMP THREADPRIVATE( icount,zday_last) 154 155 REAL zps_av 156 157 real,allocatable,SAVE :: tsurf(:),tsoil(:,:),rnatur(:) 158 real,allocatable,SAVE :: capcal(:),fluxgrd(:) 159 real,allocatable,SAVE :: dtrad(:,:),fluxrad(:) 160 real,allocatable,SAVE :: q2(:,:),q2l(:,:) 161 real,allocatable,SAVE :: albedo(:),emissiv(:),z0(:),inertie(:) 162 real,SAVE :: solarcst=1370. 163 real,SAVE :: stephan=5.67e-08 164 165 !$OMP THREADPRIVATE(tsurf,tsoil,rnatur) 166 !$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad) 167 !$OMP THREADPRIVATE( q2,q2l) 168 !$OMP THREADPRIVATE( albedo,emissiv,solarcst,z0,inertie) 169 !$OMP THREADPRIVATE( stephan) 170 171 172 EXTERNAL ismin,ismax 173 174 175 INTEGER longcles 176 PARAMETER ( longcles = 20 ) 177 REAL clesphy0( longcles ) 178 REAL presnivs(nlayer) 179 180 print*,'OK DANS PHYPARAM' 181 182 !----------------------------------------------------------------------- 183 ! 1. Initialisations : 184 ! -------------------- 185 186 ! call initial0(ngrid*nlayermx*nqmx,pdq) 187 nlevel=nlayer+1 188 189 ! print*,'OK ',nlevel 190 191 igout=ngrid/2+1 192 ! print*,'OK PHYPARAM ',ngrid,igout 193 194 195 zday=rjourvrai+gmtime 196 197 ! print*,'OK PHYPARAM 0A nq ',nq 198 tauq(1)=1800. 199 tauq(2)=10800. 200 tauq(3)=86400. 201 tauq(4)=864000. 202 ! print*,'OK PHYPARAM 0 B nq ',nq 203 factq(1:4)=(1.-exp(-ptimestep/tauq(1:4)))/ptimestep 204 205 ! print*,'OK PHYPARAM 0 ' 206 print*,'nq ',nq 207 print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid) 208 print*,'nlayer',nlayer 209 print*,'size pdq ',ngrid*nlayer*4,ngrid*nlayer*nq, & 210 & size(pdq),size(lati),size(pq),size(factq) 211 do iq=1,4 212 do l=1,nlayer 213 pdq(1:ngrid,l,iq)= & 214 & (1.+sin(lati(1:ngrid))-pq(1:ngrid,l,iq))*factq(iq) 215 enddo 216 enddo 217 218 IF(firstcall) THEN 219 220 print*,'AKk',ngrid,nsoilmx 221 allocate(tsurf(ngrid)) 222 print*,'AKa' 223 allocate (tsoil(ngrid,nsoilmx)) 224 print*,'AKb' 225 allocate (rnatur(ngrid)) 226 print*,'AK2' 227 allocate(capcal(ngrid),fluxgrd(ngrid)) 228 print*,'AK3' 229 allocate(dtrad(ngrid,nlayer),fluxrad(ngrid)) 230 print*,'AK4' 231 allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1)) 232 print*,'AK5' 233 allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid)) 234 print*,'AK6' 235 236 237 do l=1,nlayer 238 pdq(:,l,5)=1.+sin(lati(:))/ptimestep 239 enddo 240 PRINT*,'FIRSTCALL ' 241 242 ! zday_last=rjourvrai 243 zday_last=zday-ptimestep/unjours 244 ! CALL readfi(ngrid,nlayer,nsoilmx,ldrs, 245 ! . day_ini,gmtime,albedo,inertie,emissiv,z0,rnatur, 246 ! . q2,q2l,tsurf,tsoil) 247 rnatur=1. 248 emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter 249 inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter 250 albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter 251 z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter 252 q2=1.e-10 253 q2l=1.e-10 254 tsurf=300. 255 tsoil=300. 256 print*,tsoil(ngrid/2+1,nsoilmx/2+2) 257 print*,'TS ',tsurf(igout),tsoil(igout,5) 258 CALL iniorbit 259 260 if (.not.callrad) then 261 do ig=1,ngrid 262 fluxrad(ig)=0. 263 enddo 264 endif 265 266 ! print*,'OK PHYPARAM 1 ' 267 IF(callsoil) THEN 268 CALL soil(ngrid,nsoilmx,firstcall,inertie, & 269 & ptimestep,tsurf,tsoil,capcal,fluxgrd) 270 ELSE 271 PRINT*,'WARNING!!! Thermal conduction in the soil turned off' 272 DO ig=1,ngrid 273 capcal(ig)=1.e5 274 fluxgrd(ig)=0. 275 ENDDO 276 ENDIF 277 ! CALL inifrict(ptimestep) 278 icount=0 279 280 PRINT*,'FIRSTCALL B ' 281 print*,'INIIO avant iophys_ini ' 282 call iophys_ini('phys.nc ',presnivs) 283 284 ENDIF 285 icount=icount+1 286 287 PRINT*,'FIRSTCALL AP ' 288 IF (ngrid.NE.ngridmax) THEN 289 PRINT*,'STOP in inifis' 290 PRINT*,'Probleme de dimenesions :' 291 PRINT*,'ngrid = ',ngrid 292 PRINT*,'ngridmax = ',ngridmax 293 STOP 294 ENDIF 295 296 DO l=1,nlayer 297 DO ig=1,ngrid 298 pdv(ig,l)=0. 299 pdu(ig,l)=0. 300 pdt(ig,l)=0. 301 ENDDO 302 ENDDO 303 304 DO ig=1,ngrid 305 pdpsrf(ig)=0. 306 zflubid(ig)=0. 307 zdtsrf(ig)=0. 308 ENDDO 309 310 !----------------------------------------------------------------------- 311 ! calcul du geopotentiel aux niveaux intercouches 312 ! ponderation des altitudes au niveau des couches en dp/p 313 314 DO l=1,nlayer 315 DO ig=1,ngrid 316 zzlay(ig,l)=pphi(ig,l)/g 317 ENDDO 318 ENDDO 319 DO ig=1,ngrid 320 zzlev(ig,1)=0. 321 ENDDO 322 DO l=2,nlayer 323 DO ig=1,ngrid 324 z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) 325 z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) 326 zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) 327 ENDDO 328 ENDDO 329 330 !----------------------------------------------------------------------- 331 ! Transformation de la temperature en temperature potentielle 332 DO l=1,nlayer 333 DO ig=1,ngrid 334 zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp 335 ENDDO 336 ENDDO 337 DO l=1,nlayer 338 DO ig=1,ngrid 339 zh(ig,l)=pt(ig,l)/zpopsk(ig,l) 340 ENDDO 341 ENDDO 342 343 !----------------------------------------------------------------------- 344 ! 2. Calcul of the radiative tendencies : 345 ! --------------------------------------- 346 347 ! print*,'callrad0' 348 IF(callrad) THEN 349 ! print*,'callrad' 350 351 ! WARNING !!! on calcule le ray a chaque appel 352 ! IF( MOD(icount,iradia).EQ.0) THEN 353 354 CALL solarlong(zday,zls) 355 CALL orbite(zls,dist_sol,declin) 356 ! declin=0. 357 ! print*,'ATTENTIOn : pdeclin = 0',' L_s=',zls 358 print*,'diurnal=',diurnal 359 IF(diurnal) THEN 360 if ( 1.eq.1 ) then 361 ztim1=SIN(declin) 362 ztim2=COS(declin)*COS(2.*pi*(zday-.5)) 363 ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) 364 ! call dump2d(iim,jjm-1,lati(2),'LATI 0 ') 365 ! call dump2d(iim,jjm-1,long(2),'LONG 0 ') 366 ! call dump2d(iim,jjm-1,sinlon(2),'sinlon0 ') 367 ! call dump2d(iim,jjm-1,coslon(2),'coslon0 ') 368 ! call dump2d(iim,jjm-1,sinlat(2),'sinlat ') 369 ! call dump2d(iim,jjm-1,coslat(2),'coslat ') 370 CALL solang(ngrid,sinlon,coslon,sinlat,coslat, & 371 & ztim1,ztim2,ztim3, & 372 & mu0,fract) 373 else 374 zdtime=ptimestep*float(iradia) 375 CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract) 376 print*,'ZENANG ' 377 endif 378 379 IF(lverbose) THEN 380 PRINT*,'day, declin, sinlon,coslon,sinlat,coslat' 381 PRINT*,zday, declin, & 382 & sinlon(igout),coslon(igout), & 383 & sinlat(igout),coslat(igout) 384 ENDIF 385 ELSE 386 print*,'declin,ngrid,rad',declin,ngrid,rad 387 388 ! call dump2d(iim,jjm-1,lati(2),'LATI ') 389 CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad) 390 ENDIF 391 ! call dump2d(iim,jjm-1,fract(2),'FRACT A ') 392 ! call dump2d(iim,jjm-1,mu0(2),'MU0 A ') 393 394 395 ! 2.2 Calcul of the radiative tendencies and fluxes: 396 ! -------------------------------------------------- 397 398 ! 2.1.2 levels 399 400 zinsol=solarcst/(dist_sol*dist_sol) 401 print*,iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer 402 print*,'iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer' 403 ! call dump2d(iim,jjm-1,albedo(2),'ALBEDO ') 404 ! call dump2d(iim,jjm-1,mu0(2),'MU0 ') 405 ! call dump2d(iim,jjm-1,fract(2),'FRACT ') 406 ! call dump2d(iim,jjm-1,lati(2),'LATI ') 407 zps_av=1.e5 408 if (firstcall) then 409 print*,'WARNING ps_rad impose' 410 endif 411 CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, & 412 & pplev,zps_av, & 413 & mu0,fract,zinsol, & 414 & zfluxsw,zdtsw, & 415 & lverbose) 416 ! call dump2d(iim,jjm-1,zfluxsw(2),'SWS 1 ') 417 ! stop 418 419 CALL lw(ngrid,nlayer,coefir,emissiv, & 420 & pplev,zps_av,tsurf,pt, & 421 & zfluxlw,zdtlw, & 422 & lverbose) 423 424 425 ! 2.4 total flux and tendencies: 426 ! ------------------------------ 427 428 ! 2.4.1 fluxes 429 430 DO ig=1,ngrid 431 fluxrad(ig)=emissiv(ig)*zfluxlw(ig) & 432 & +zfluxsw(ig)*(1.-albedo(ig)) 433 zplanck(ig)=tsurf(ig)*tsurf(ig) 434 zplanck(ig)=emissiv(ig)* & 435 & stephan*zplanck(ig)*zplanck(ig) 436 fluxrad(ig)=fluxrad(ig)-zplanck(ig) 437 ENDDO 438 439 ! 2.4.2 temperature tendencies 440 441 DO l=1,nlayer 442 DO ig=1,ngrid 443 dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) 444 ENDDO 445 ENDDO 446 447 ! ENDIF 448 449 ! 2.5 Transformation of the radiative tendencies: 450 ! ----------------------------------------------- 451 452 DO l=1,nlayer 453 DO ig=1,ngrid 454 pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) 455 ENDDO 456 ENDDO 457 458 IF(lverbose) THEN 459 PRINT*,'Diagnotique for the radaition' 460 PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck' 461 PRINT*,albedo(igout),emissiv(igout),mu0(igout), & 462 & fract(igout), & 463 & fluxrad(igout),zplanck(igout) 464 PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)' 465 PRINT*,'unjours',unjours 466 DO l=1,nlayer 467 WRITE(*,'(3f15.5,2e15.2)') pt(igout,l), & 468 & pplay(igout,l),pplev(igout,l), & 469 & zdtsw(igout,l),zdtlw(igout,l) 470 ENDDO 471 ENDIF 472 473 ENDIF 474 475 !----------------------------------------------------------------------- 476 ! 3. Vertical diffusion (turbulent mixing): 477 ! ----------------------------------------- 478 ! 479 IF(calldifv) THEN 480 481 DO ig=1,ngrid 482 zflubid(ig)=fluxrad(ig)+fluxgrd(ig) 483 ENDDO 484 485 zdum1(:,:)=0. 486 zdum2(:,:)=0. 487 488 do l=1,nlayer 489 do ig=1,ngrid 490 zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l) 491 enddo 492 enddo 493 494 CALL vdif(ngrid,nlayer,zday, & 495 & ptimestep,capcal,z0, & 496 & pplay,pplev,zzlay,zzlev, & 497 & pu,pv,zh,tsurf,emissiv, & 498 & zdum1,zdum2,zdum3,zflubid, & 499 & zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l, & 500 & lverbose) 501 502 DO l=1,nlayer 503 DO ig=1,ngrid 504 pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l) 505 pdu(ig,l)=pdu(ig,l)+zdufr(ig,l) 506 pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l) 507 ENDDO 508 ENDDO 509 510 DO ig=1,ngrid 511 zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig) 512 ENDDO 513 514 ELSE 515 DO ig=1,ngrid 516 zdtsrf(ig)=zdtsrf(ig)+ & 517 & (fluxrad(ig)+fluxgrd(ig))/capcal(ig) 518 ENDDO 519 ENDIF 520 ! 521 !----------------------------------------------------------------------- 522 ! 4. Dry convective adjustment: 523 ! ----------------------------- 524 525 IF(calladj) THEN 526 527 DO l=1,nlayer 528 DO ig=1,ngrid 529 zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l) 530 ENDDO 531 ENDDO 532 533 zdufr(:,:)=0. 534 zdvfr(:,:)=0. 535 zdhfr(:,:)=0. 536 537 CALL convadj(ngrid,nlayer,ptimestep, & 538 & pplay,pplev,zpopsk, & 539 & pu,pv,zh, & 540 & pdu,pdv,zdum1, & 541 & zdufr,zdvfr,zdhfr) 542 543 DO l=1,nlayer 544 DO ig=1,ngrid 545 pdu(ig,l)=pdu(ig,l)+zdufr(ig,l) 546 pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l) 547 pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l) 548 ENDDO 549 ENDDO 550 551 ENDIF 552 553 !----------------------------------------------------------------------- 554 ! On incremente les tendances physiques de la temperature du sol: 555 ! --------------------------------------------------------------- 556 557 DO ig=1,ngrid 558 tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) 559 ENDDO 560 561 WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1) 562 563 !----------------------------------------------------------------------- 564 ! soil temperatures: 565 ! -------------------- 566 567 IF (callsoil) THEN 568 CALL soil(ngrid,nsoilmx,.false.,inertie, & 569 & ptimestep,tsurf,tsoil,capcal,fluxgrd) 570 IF(lverbose) THEN 571 PRINT*,'Surface Heat capacity,conduction Flux, Ts, dTs, dt' 572 PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout), & 573 & zdtsrf(igout),ptimestep 574 ENDIF 575 ENDIF 576 577 !----------------------------------------------------------------------- 578 ! sorties: 579 ! -------- 580 581 ! call dump2d(iim,jjm-1,zfluxsw(2),'SWS 2 ') 582 print*,'zday, zday_last ',zday,zday_last,icount 583 if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then 584 print*,'zday, zday_last SORTIE ',zday,zday_last 604 585 zday_last=zday 605 cEcriture/extension de la coordonnee temps606 607 608 609 610 586 ! Ecriture/extension de la coordonnee temps 587 588 do ig=1,ngridmax 589 zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*285.)) 590 enddo 591 611 592 call iophys_ecrit('u',llm,'Vent zonal moy','m/s',pu) 612 593 call iophys_ecrit('v',llm,'Vent meridien moy','m/s',pv) … … 614 595 call iophys_ecrit('geop',llm,'Geopotential','m2/s2',pphi) 615 596 call iophys_ecrit('plev',llm,'plev','Pa',pplev(:,1:nlayer)) 616 597 617 598 call iophys_ecrit('du',llm,'du',' ',pdu) 618 599 call iophys_ecrit('dv',llm,'du',' ',pdv) … … 620 601 call iophys_ecrit('dtsw',llm,'dtsw',' ',zdtsw) 621 602 call iophys_ecrit('dtlw',llm,'dtlw',' ',zdtlw) 622 603 623 604 do iq=1,nq 624 605 nomq="tr." … … 626 607 call iophys_ecrit(nomq,llm,nomq,' ',pq(:,:,iq)) 627 608 enddo 628 609 629 610 call iophys_ecrit('ts',1,'Surface temper','K',tsurf) 630 611 call iophys_ecrit('coslon',1,'coslon',' ',coslon) … … 639 620 call iophys_ecrit('swsurf',1,'SW surf','Pa',zfluxsw) 640 621 call iophys_ecrit('lwsurf',1,'LW surf','Pa',zfluxlw) 641 642 643 644 c-----------------------------------------------------------------------645 646 647 648 ! if (ierr.ne.0) then649 ! write(6,*)' Pb d''ouverture du fichier restart'650 ! write(6,*)' ierr = ', ierr651 ! call exit(1)652 ! endif653 654 cs (tsurf(1),ig=1,iim+1),655 cs ( (tsurf(ig),ig=(j-2)*iim+2,(j-1)*iim+1),656 cs tsurf((j-2)*iim+2) ,j=2,jjm),657 cs (tsurf(ngridmax),ig=1,iim+1),658 cs ( (tsoil(1,l),ig=1,iim+1),659 cs ( (tsoil(ig,l),ig=(j-2)*iim+2,(j-1)*iim+1),660 cs tsoil((j-2)*iim+2,l) ,ig=2,jjm),661 cs (tsoil(ngridmax,l),ig=1,iim+1)662 cs ,l=1,nsoilmx)663 ENDIF664 665 666 RETURN 667 END 622 623 endif 624 625 !----------------------------------------------------------------------- 626 IF(lastcall) THEN 627 call iotd_fin 628 PRINT*,'Ecriture du fichier de reinitialiastion de la physique' 629 ! if (ierr.ne.0) then 630 ! write(6,*)' Pb d''ouverture du fichier restart' 631 ! write(6,*)' ierr = ', ierr 632 ! call exit(1) 633 ! endif 634 write(75) tsurf,tsoil 635 ! s (tsurf(1),ig=1,iim+1), 636 ! s ( (tsurf(ig),ig=(j-2)*iim+2,(j-1)*iim+1), 637 ! s tsurf((j-2)*iim+2) ,j=2,jjm), 638 ! s (tsurf(ngridmax),ig=1,iim+1), 639 ! s ( (tsoil(1,l),ig=1,iim+1), 640 ! s ( (tsoil(ig,l),ig=(j-2)*iim+2,(j-1)*iim+1), 641 ! s tsoil((j-2)*iim+2,l) ,ig=2,jjm), 642 ! s (tsoil(ngridmax,l),ig=1,iim+1) 643 ! s ,l=1,nsoilmx) 644 ENDIF 645 646 END SUBROUTINE phyparam 647 648 END MODULE phyparam_mod -
dynamico_lmdz/simple_physics/phyparam/param/physiq_mod.F90
r4211 r4212 28 28 29 29 USE infotrac_phy, only : nqtot 30 USE phyparam_mod, ONLY : phyparam 30 31 ! 31 32 ! Routine argument: … … 52 53 real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure 53 54 54 real :: temp_newton(klon,klev)55 integer :: k56 55 logical, save :: first=.true. 57 56 !$OMP THREADPRIVATE(first) … … 61 60 REAL,SAVE :: gmtime=0. 62 61 !$OMP THREADPRIVATE(rjourvrai,gmtime) 63 64 62 65 63 ! initializations 66 64 if (debut) CALL init_physiq_mod(pdtphys, presnivs) ! Things to do only for the first call to physics
Note: See TracChangeset
for help on using the changeset viewer.
