![(please configure the [header_logo] section in trac.ini)](/dynamico/chrome/site/your_project_logo.png){kind=logo}
Changeset 186 for codes/icosagcm/trunk/src/phyparam.F
- Timestamp:
- 01/09/14 09:56:11 (10 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/phyparam.F (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/phyparam.F
r163 r186 1 MODULE PHY2 USE dimphys_mod3 LOGICAL:: firstcall,lastcall4 contains5 6 SUBROUTINE phyparam_lmd(it,ngrid,nlayer,nq,7 s ptimestep,lati,long,rjourvrai,gmtime,8 s pplev,pplay,pphi,pphis,9 s pu,pv,pt,pq,10 s pdu,pdv,pdt,pdq,pdpsrf)11 12 USE ICOSA13 USE dimphys_mod14 USE RADIATION15 USE SURFACE_PROCESS16 c17 IMPLICIT NONE18 c=======================================================================19 c20 c subject:21 c --------22 c23 c Organisation of the physical parametrisations of the LMD24 c 20 parameters GCM for planetary atmospheres.25 c It includes:26 c raditive transfer (long and shortwave) for CO2 and dust.27 c vertical turbulent mixing28 c convective adjsutment29 c30 c author: Frederic Hourdin 15 / 10 /9331 c -------32 c33 c arguments:34 c ----------35 c36 c input:37 c ------38 c39 c ngrid Size of the horizontal grid.40 c All internal loops are performed on that grid.41 c nlayer Number of vertical layers.42 c nq Number of advected fields43 c firstcall True at the first call44 c lastcall True at the last call45 c rjourvrai Number of days counted from the North. Spring46 c equinoxe.47 c gmtime hour (s)48 c ptimestep timestep (s)49 c pplay(ngrid,nlayer+1) Pressure at the middle of the layers (Pa)50 c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)51 c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)52 c pu(ngrid,nlayer) u component of the wind (ms-1)53 c pv(ngrid,nlayer) v component of the wind (ms-1)54 c pt(ngrid,nlayer) Temperature (K)55 c pq(ngrid,nlayer,nq) Advected fields56 c pudyn(ngrid,nlayer) \57 c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the58 c ptdyn(ngrid,nlayer) / corresponding variables59 c pqdyn(ngrid,nlayer,nq) /60 c pw(ngrid,?) vertical velocity61 c62 c output:63 c -------64 c65 c pdu(ngrid,nlayer) \66 c pdv(ngrid,nlayer) \ Temporal derivative of the corresponding67 c pdt(ngrid,nlayer) / variables due to physical processes.68 c pdq(ngrid,nlayer) /69 c pdpsrf(ngrid) /70 c71 c=======================================================================72 c73 !c-----------------------------------------------------------------------74 !c75 !c 0. Declarations :76 !c ------------------77 !c Arguments :78 !c -----------79 80 !c inputs:81 !c -------82 INTEGER ngrid,nlayer,nq,it,ij,i83 REAL ptimestep84 real zdtime85 REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)86 REAL pphi(ngrid,nlayer)87 REAL pphis(ngrid)88 REAL pu(ngrid,nlayer),pv(ngrid,nlayer)89 REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)90 REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)91 92 !c dynamial tendencies93 REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq)94 REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer)95 REAL pw(ngrid,nlayer)96 97 !c Time98 real rjourvrai99 REAL gmtime100 101 !c outputs:102 !c --------103 104 !c physical tendencies105 REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)106 REAL pdpsrf(ngrid)107 108 109 !c Local variables :110 !c -----------------111 112 INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,unit,isoil113 REAL zps_av114 REAL zday115 REAL zh(ngrid,nlayer),z1,z2116 REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)117 REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)118 REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)119 REAL zflubid(ngrid),zpmer(ngrid)120 REAL zplanck(ngrid),zpopsk(ngrid,nlayer)121 REAL zdum1(ngrid,nlayer)122 REAL zdum2(ngrid,nlayer)123 REAL zdum3(ngrid,nlayer)124 REAL ztim1,ztim2,ztim3125 REAL zls,zinsol126 REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)127 REAL zfluxsw(ngrid),zfluxlw(ngrid)128 character*2 str2129 130 !c Local saved variables:131 !c ----------------------132 REAL(rstd)::long(ngrid),lati(ngrid)133 REAL(rstd)::mu0(ngrid),fract(ngrid),coslat(ngrid)134 REAL(rstd)::sinlon(ngrid),coslon(ngrid),sinlat(ngrid)135 REAL(rstd)::dist_sol,declin136 REAL::totarea !sarvesh137 !!!!!!!!sarvesh !!!!!!! CHECK SAVE ATTRIBUTE138 INTEGER:: icount139 real:: zday_last140 REAL:: solarcst141 142 SAVE icount,zday_last143 SAVE solarcst144 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!145 REAL stephan146 SAVE stephan147 DATA stephan/5.67e-08/148 DATA solarcst/1370./149 REAL presnivs(nlayer)150 INTEGER:: nn1,nn2151 c-----------------------------------------------------------------------152 c 1. Initialisations :153 c --------------------154 c call initial0(ngrid*nlayer*nqmx,pdq)155 156 ! nn1=(jj_begin -1)*iim+ii_begin157 ! nn2=(jj_end -1)*iim+ii_end158 159 nlevel=nlayer+1160 igout=ngrid/2+1161 162 DO j=jj_begin-offset,jj_end+offset163 DO i=ii_begin-offset,ii_end+offset164 ig=(j-1)*iim+i165 sinlat(ig) = sin(lati(ig))166 coslat(ig) = cos(lati(ig))167 sinlon(ig) = sin(long(ig))168 coslon(ig) = cos(long(ig))169 END DO170 ENDDO171 zday=rjourvrai+gmtime172 IF ( it == 0 ) then173 firstcall=.TRUE.174 ELSE175 firstcall=.FALSE.176 ENDIF177 IF ( it == ndays*day_step ) Then178 lastcall = .True.179 END IF180 181 IF(firstcall) THEN182 PRINT*,'FIRSTCALL ',ngridmx,nlayermx,nsoilmx183 zday_last=rjourvrai184 inertie=2000185 albedo=0.2186 emissiv=1.187 z0=0.1188 rnatur=1.189 q2=1.e-10190 q2l=1.e-10191 tsurf(:)=300.192 tsoil(:,:)=300.193 ! print*,tsoil(ngrid/2+1,nsoilmx/2+2)194 ! print*,'TS ',tsurf(igout),tsoil(igout,5)195 196 IF (.not.callrad) then197 DO j=jj_begin-offset,jj_end+offset198 DO i=ii_begin-offset,ii_end+offset199 ig=(j-1)*iim+i200 fluxrad(ig)=0.201 enddo202 enddo203 ENDIF204 205 IF(callsoil) THEN206 CALL soil(ngrid,nsoilmx,firstcall,inertie,207 s ptimestep,tsurf,tsoil,capcal,fluxgrd)208 ELSE209 PRINT*,'WARNING!!! Thermal conduction in the soil210 s turned off'211 DO j=jj_begin-offset,jj_end+offset212 DO i=ii_begin-offset,ii_end+offset213 ig=(j-1)*iim+i214 capcal(ig)=1.e5215 fluxgrd(ig)=0.216 ENDDO217 ENDDO218 ENDIF219 icount=0220 ENDIF221 222 IF (ngrid.NE.ngrid) THEN223 PRINT*,'STOP in inifis'224 PRINT*,'Probleme de dimenesions :'225 PRINT*,'ngrid = ',ngrid226 PRINT*,'ngrid = ',ngrid227 STOP228 ENDIF229 230 DO l=1,nlayer231 DO j=jj_begin-offset,jj_end+offset232 DO i=ii_begin-offset,ii_end+offset233 ig=(j-1)*iim+i234 pdv(ig,l)=0.235 pdu(ig,l)=0.236 pdt(ig,l)=0.237 ENDDO238 ENDDO239 ENDDO240 241 DO j=jj_begin-offset,jj_end+offset242 DO i=ii_begin-offset,ii_end+offset243 ig=(j-1)*iim+i244 pdpsrf(ig)=0.245 zflubid(ig)=0.246 zdtsrf(ig)=0.247 ENDDO248 ENDDO249 250 zps_av=0.0251 DO j=jj_begin-offset,jj_end+offset252 DO i=ii_begin-offset,ii_end+offset253 ig=(j-1)*iim+i254 zps_av=zps_av+pplev(ig,1)*Ai(ig)255 totarea=totarea+Ai(ig)256 END DO257 END DO258 zps_av=zps_av/totarea259 260 261 !print*,"maxplev",maxval(pplev(:,1)),minval(pplev(:,1))262 c-----------------------------------------------------------------------263 c calcul du geopotentiel aux niveaux intercouches264 c ponderation des altitudes au niveau des couches en dp/p265 266 DO l=1,nlayer267 DO j=jj_begin-offset,jj_end+offset268 DO i=ii_begin-offset,ii_end+offset269 ig=(j-1)*iim+i270 zzlay(ig,l)=pphi(ig,l)/g271 ENDDO272 ENDDO273 ENDDO274 275 !print*,"zzlay",maxval(zzlay(:,1)),minval(zzlay(:,1))276 277 278 DO j=jj_begin-offset,jj_end+offset279 DO i=ii_begin-offset,ii_end+offset280 ig=(j-1)*iim+i281 zzlev(ig,1)=0.282 ENDDO283 ENDDO284 285 DO l=2,nlayer286 DO j=jj_begin-offset,jj_end+offset287 DO i=ii_begin-offset,ii_end+offset288 ig=(j-1)*iim+i289 z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))290 z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))291 zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)292 ENDDO293 ENDDO294 ENDDO295 !print*,"zzlev",maxval(zzlev(:,1)),minval(zzlev(:,1))296 c-----------------------------------------------------------------------297 c Transformation de la temperature en temperature potentielle298 DO l=1,nlayer299 DO j=jj_begin-offset,jj_end+offset300 DO i=ii_begin-offset,ii_end+offset301 ig=(j-1)*iim +i302 zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**kappa303 ENDDO304 ENDDO305 ENDDO306 307 DO l=1,nlayer308 DO j=jj_begin-offset,jj_end+offset309 DO i=ii_begin-offset,ii_end+offset310 ig=(j-1)*iim+i311 zh(ig,l)=pt(ig,l)/zpopsk(ig,l)312 ENDDO313 ENDDO314 ENDDO315 !print*,"ph pot",maxval(zh(:,1)),minval(zh(:,1))316 ! go to 101317 c-----------------------------------------------------------------------318 c 2. Calcul of the radiative tendencies :319 c ---------------------------------------320 ! print*,'callrad0'321 IF(callrad) THEN322 ! print*,'callrad'323 ! WARNING !!! on calcule le ray a chaque appel324 ! IF( MOD(icount,iradia).EQ.0) THEN325 326 CALL solarlong(zday,zls)327 CALL orbite(zls,dist_sol,declin)328 !329 ! print*,'ATTENTIOn : pdeclin = 0',' L_s=',zls330 ! print*,'diurnal=',diurnal331 IF(diurnal) THEN332 if ( 1.eq.1 ) then333 ztim1=SIN(declin)334 ztim2=COS(declin)*COS(2.*pi*(zday-.5))335 ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))336 337 CALL solang(ngrid,sinlon,coslon,sinlat,coslat,338 s ztim1,ztim2,ztim3,339 s mu0,fract)340 else341 zdtime=ptimestep*float(iradia)342 ! CALL zenang(zls,gmtime,zdtime,lati,long,mu0,fract) ! FIXME343 !print*,'ZENANG '344 endif345 346 IF(lverbose) THEN347 PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'348 PRINT*,zday, declin,349 s sinlon(igout),coslon(igout),350 s sinlat(igout),coslat(igout)351 ENDIF352 ELSE353 !print*,'declin,ngrid,radius',declin,ngrid,radius354 CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,radius)355 ! open(100,file="mu0.txt")356 ! write(100,*)(mu0(ij),ij=1,ngrid)357 ENDIF358 ! print*,"iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii"359 !c 2.2 Calcul of the radiative tendencies and fluxes:360 !c --------------------------------------------------361 !c 2.1.2 levels362 363 zinsol=solarcst/(dist_sol*dist_sol)364 ! print*,'iim,jjm,llm,ngrid,nlayer,ngrid,nlayer'365 ! print*,iim,jjm,llm,ngrid,nlayer,ngrid,nlayer366 ! print*,"zinsol sol_dist",zinsol,dist_sol367 ! STOP368 369 CALL sw(ngrid,nlayer,diurnal,coefvis,albedo,370 $ pplev,zps_av,371 $ mu0,fract,zinsol,372 $ zfluxsw,zdtsw,373 $ lverbose)374 375 !print*,"sw",maxval(zfluxsw),minval(zfluxsw),376 ! $ maxval(zdtsw),minval(zdtsw), " it",it377 378 ! print*,"lllllllllllllllllllllllllllllllllllllllll"379 ! print*,"pplev",maxval(pplev),minval(pplev)380 ! print*,"zps,tsurf",zps_av,maxval(tsurf),minval(tsurf)381 ! print*,"pt",maxval(pt),minval(pt)382 383 384 CALL lw(ngrid,nlayer,coefir,emissiv,385 $ pplev,zps_av,tsurf,pt,386 $ zfluxlw,zdtlw,387 $ lverbose)388 389 !print*,"lw",maxval(zfluxlw),minval(zfluxlw),390 ! $ maxval(zdtlw),minval(zdtlw)391 392 ! print*,"lw",maxval(393 394 ! print*,"after lw"395 c 2.4 total flux and tendencies:396 c ------------------------------397 398 c 2.4.1 fluxes399 400 DO j=jj_begin-offset,jj_end+offset401 DO i=ii_begin-offset,ii_end+offset402 ig=(j-1)*iim+i403 fluxrad(ig)=emissiv(ig)*zfluxlw(ig)404 $ +zfluxsw(ig)*(1.-albedo(ig))405 406 zplanck(ig)=tsurf(ig)*tsurf(ig)407 408 zplanck(ig)=emissiv(ig)*409 $ stephan*zplanck(ig)*zplanck(ig)410 411 fluxrad(ig)=fluxrad(ig)-zplanck(ig)412 ENDDO413 ENDDO414 c 2.4.2 temperature tendencies415 416 DO l=1,nlayer417 DO j=jj_begin-offset,jj_end+offset418 DO i=ii_begin-offset,ii_end+offset419 ig=(j-1)*iim+i420 dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)421 ENDDO422 ENDDO423 ENDDO424 425 426 c 2.5 Transformation of the radiative tendencies:427 c -----------------------------------------------428 429 DO l=1,nlayer430 DO j=jj_begin-offset,jj_end+offset431 DO i=ii_begin-offset,ii_end+offset432 ig=(j-1)*iim+i433 pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)434 ENDDO435 ENDDO436 ENDDO437 438 IF(lverbose) THEN439 PRINT*,'Diagnotique for the radaition'440 PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'441 PRINT*,albedo(igout),emissiv(igout),mu0(igout),442 s fract(igout),443 s fluxrad(igout),zplanck(igout)444 PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'445 ! PRINT*,'unjours',unjours446 DO l=1,nlayer447 WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),448 s pplay(igout,l),pplev(igout,l),449 s zdtsw(igout,l),zdtlw(igout,l)450 ENDDO451 ENDIF452 ENDIF !( CALL RADIATION )453 ! print*,"eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee"454 455 c-----------------------------------------------------------------------456 c 3. Vertical diffusion (turbulent mixing):457 c -----------------------------------------458 c459 IF(calldifv) THEN460 461 DO ig=1,ngrid462 zflubid(ig)=fluxrad(ig)+fluxgrd(ig)463 ENDDO464 465 CALL zerophys(ngrid*nlayer,zdum1)466 CALL zerophys(ngrid*nlayer,zdum2)467 c CALL zerophys(ngrid*nlayer,zdum3)468 do l=1,nlayer469 do ig=1,ngrid470 zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l)471 enddo472 enddo473 474 CALL vdif(ngrid,nlayer,zday,475 $ ptimestep,capcal,z0,476 $ pplay,pplev,zzlay,zzlev,477 $ pu,pv,zh,tsurf,emissiv,478 $ zdum1,zdum2,zdum3,zflubid,479 $ zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l,480 $ lverbose)481 c igout=ngrid/2+1482 c PRINT*,'zdufr zdvfr zdhfr'483 c DO l=1,nlayer484 c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)485 c ENDDO486 c CALL difv (ngrid,nlayer,487 c $ area,lati,capcal,488 c $ pplev,pphi,489 c $ pu,pv,zh,tsurf,emissiv,490 c $ zdum1,zdum2,zdum3,zflubid,491 c $ zdufr,zdvfr,zdhfr,zdtsrf,492 c $ .true.)493 c PRINT*,'zdufr zdvfr zdhfr'494 c DO l=1,nlayer495 c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l)496 c ENDDO497 c STOP498 499 DO l=1,nlayer500 DO ig=1,ngrid501 pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)502 pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)503 pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)504 ENDDO505 ENDDO506 507 DO ig=1,ngrid508 zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig)509 ENDDO510 511 ELSE512 513 DO j=jj_begin-offset,jj_end+offset514 DO i=ii_begin-offset,ii_end+offset515 ig=(j-1)*iim+i516 zdtsrf(ig)=zdtsrf(ig)+517 s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)518 ENDDO519 ENDDO520 521 ! write(66,*)"tsrf",maxval(zdtsrf(:)),minval(zdtsrf(:))522 ! write(66,*)"frd",maxval(fluxrad(:)),minval(fluxrad(:))523 ! write(66,*)"fgd",maxval(fluxgrd(:)),minval(fluxgrd(:))524 ENDIF525 c-----------------------------------------------------------------------526 c 4. Dry convective adjustment:527 c -----------------------------528 529 IF(calladj) THEN530 531 DO l=1,nlayer532 DO ig=1,ngrid533 zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l)534 ENDDO535 ENDDO536 CALL zerophys(ngrid*nlayer,zdufr)537 CALL zerophys(ngrid*nlayer,zdvfr)538 CALL zerophys(ngrid*nlayer,zdhfr)539 CALL convadj(ngrid,nlayer,ptimestep,540 $ pplay,pplev,zpopsk,541 $ pu,pv,zh,542 $ pdu,pdv,zdum1,543 $ zdufr,zdvfr,zdhfr)544 545 DO l=1,nlayer546 DO ig=1,ngrid547 pdu(ig,l)=pdu(ig,l)+zdufr(ig,l)548 pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l)549 pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l)550 ENDDO551 ENDDO552 553 ENDIF554 555 !101 continue556 c-----------------------------------------------------------------------557 c On incremente les tendances physiques de la temperature du sol:558 c ---------------------------------------------------------------559 ! WRITE(55,*)"tsurf",maxval(tsurf(:)),minval(tsurf(:)),it560 561 DO j=jj_begin-offset,jj_end+offset562 DO i=ii_begin-offset,ii_end+offset563 ig=(j-1)*iim+i564 tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)565 ENDDO566 ENDDO567 c-----------------------------------------------------------------------568 c soil temperatures:569 c --------------------570 571 IF (callsoil) THEN572 CALL soil(ngrid,nsoilmx,.false.,inertie,573 s ptimestep,tsurf,tsoil,capcal,fluxgrd)574 IF(lverbose) THEN575 PRINT*,'Surface Heat capacity,conduction Flux, Ts,576 s dTs, dt'577 PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout),578 s zdtsrf(igout),ptimestep579 ENDIF580 ENDIF581 582 c-----------------------------------------------------------------------583 c sorties:584 c --------585 if(zday.GT.zday_last+period_sort) then586 zday_last=zday587 c Ecriture/extension de la coordonnee temps588 589 do ig=1,ngrid590 zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(kappa*cpp*285.))591 enddo592 endif593 c-----------------------------------------------------------------------594 IF(lastcall) THEN595 PRINT*,'Ecriture du fichier de reinitialiastion de la physique'596 ENDIF597 598 icount=icount+1599 ! RETURN600 END SUBROUTINE phyparam_lmd601 END MODULE PHY
Note: See TracChangeset
for help on using the changeset viewer.
