![(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/radiation_mod.F
- Timestamp:
- 01/09/14 09:56:11 (10 years ago)
- File:
-
- 1 edited
-
codes/icosagcm/trunk/src/radiation_mod.F (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/radiation_mod.F
r149 r186 1 MODULE RADIATION2 USE ICOSA3 USE dimphys_mod4 ! USE PHY5 contains6 7 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%8 SUBROUTINE zerophys(n,x)9 IMPLICIT NONE10 INTEGER n11 REAL x(n)12 INTEGER i13 14 DO i=1,n15 x(i)=0.16 ENDDO17 RETURN18 END subroutine zerophys19 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%20 21 SUBROUTINE solarlong(pday,psollong)22 IMPLICIT NONE23 c=======================================================================24 c25 c Objet:26 c ------27 c28 c Calcul de la distance soleil-planete et de la declinaison29 c en fonction du jour de l'annee.30 c31 c32 c Methode:33 c --------34 c35 c Calcul complet de l'elipse36 c37 c Interface:38 c ----------39 c40 c Uncommon comprenant les parametres orbitaux.41 c42 c Arguments:43 c ----------44 c45 c Input:46 c ------47 c pday jour de l'annee (le jour 0 correspondant a l'equinoxe)48 c lwrite clef logique pour sorties de controle49 c50 c Output:51 c -------52 c pdist_sol distance entre le soleil et la planete53 c ( en unite astronomique pour utiliser la constante54 c solaire terrestre 1370 Wm-2 )55 c pdecli declinaison ( en radians )56 c57 c=======================================================================58 c-----------------------------------------------------------------------59 c Declarations:60 c -------------61 62 c arguments:63 c ----------64 65 REAL pday,pdist_sol,pdecli,psollong66 LOGICAL lwrite67 68 c Local:69 c ------70 71 REAL zanom,xref,zx0,zdx,zteta,zz72 INTEGER iter73 74 75 c-----------------------------------------------------------------------76 c calcul de l'angle polaire et de la distance au soleil :77 c -------------------------------------------------------78 79 c calcul de l'zanomalie moyenne80 81 zz=(pday-peri_day)/year_day82 zanom=2.*pi*(zz-nint(zz))83 xref=abs(zanom)84 85 c resolution de lequation horaire zx0 - e * sin (zx0) = xref86 c methode de Newton87 88 zx0=xref+e_elips*sin(xref)89 DO 110 iter=1,1090 zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))91 if(abs(zdx).le.(1.e-7)) goto 12092 zx0=zx0+zdx93 110 continue94 120 continue95 zx0=zx0+zdx96 if(zanom.lt.0.) zx0=-zx097 98 c zteta est la longitude solaire99 100 zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))101 102 psollong=zteta-timeperi103 104 IF(psollong.LT.0.) psollong=psollong+2.*pi105 IF(psollong.GT.2.*pi) psollong=psollong-2.*pi106 107 RETURN108 END SUBROUTINE solarlong109 110 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%111 112 SUBROUTINE orbite(pls,pdist_sol,pdecli)113 IMPLICIT NONE114 c====================================================================115 c116 c Objet:117 c ------118 c119 c Distance from sun and declimation as a function of the solar120 c longitude Ls121 c122 c Interface:123 c ----------124 c Arguments:125 c ----------126 c127 c Input:128 c ------129 c pls Ls130 c131 c Output:132 c -------133 c pdist_sol Distance Sun-Planet in UA134 c pdecli declinaison ( en radians )135 c136 c====================================================================137 c-----------------------------------------------------------------------138 c Declarations:139 c -------------140 c arguments:141 c ----------142 143 REAL pday,pdist_sol,pdecli,pls144 145 c-----------------------------------------------------------------------146 147 c Distance Sun-Planet148 149 pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi))150 151 c Solar declination152 153 pdecli= asin (sin(pls)*sin(obliquit*pi/180.))154 155 c-----------------------------------------------------------------------156 c sorties eventuelles:157 c ---------------------158 159 END SUBROUTINE orbite160 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%161 162 SUBROUTINE iniorbit163 $ (paphelie,pperiheli,pyear_day,pperi_day,pobliq)164 IMPLICIT NONE165 c=====================================================================166 c167 c Auteur:168 c -------169 c Frederic Hourdin 22 Fevrier 1991170 c171 c Objet:172 c ------173 c Initialisation du sous programme orbite qui calcule174 c a une date donnee de l'annee de duree year_day commencant175 c a l'equinoxe de printemps et dont le perihelie se situe176 c a la date peri_day, la distance au soleil et la declinaison.177 c178 c Interface:179 c ----------180 c - Doit etre appele avant d'utiliser orbite.181 c - initialise le common comorbit182 c183 c Arguments:184 c ----------185 c186 c Input:187 c ------188 c aphelie \ aphelie et perihelie de l'orbite189 c periheli / en millions de kilometres.190 c191 c=====================================================================192 c Declarations:193 c -------------194 c Arguments:195 c ----------196 197 REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq198 199 c Local:200 c ------201 202 REAL zxref,zanom,zz,zx0,zdx203 INTEGER iter204 205 c'-----------------------------------------------------------------------206 207 aphelie =paphelie208 periheli=pperiheli209 year_day=pyear_day210 obliquit=pobliq211 peri_day=pperi_day212 213 PRINT*,'Perihelie en Mkm ',periheli214 PRINT*,'Aphelise en Mkm ',aphelie215 PRINT*,'obliquite en degres :',obliquit216 unitastr=149.597927217 e_elips=(aphelie-periheli)/(periheli+aphelie)218 p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr219 220 print*,'e_elips',e_elips221 print*,'p_elips',p_elips222 223 c-----------------------------------------------------------------------224 c calcul de l'angle polaire et de la distance au soleil :225 c -------------------------------------------------------226 227 c calcul de l'zanomalie moyenne228 229 zz=(year_day-pperi_day)/year_day230 zanom=2.*pi*(zz-nint(zz))231 zxref=abs(zanom)232 PRINT*,'zanom ',zanom233 234 c resolution de lequation horaire zx0 - e * sin (zx0) = zxref235 c methode de Newton236 237 zx0=zxref+e_elips*sin(zxref)238 DO 110 iter=1,100239 zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))240 if(abs(zdx).le.(1.e-12)) goto 120241 zx0=zx0+zdx242 110 continue243 120 continue244 zx0=zx0+zdx245 if(zanom.lt.0.) zx0=-zx0246 PRINT*,'zx0 ',zx0247 248 c zteta est la longitude solaire249 250 timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))251 PRINT*,'longitude solaire du perihelie timeperi = ',timeperi252 253 RETURN254 END SUBROUTINE iniorbit255 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%256 257 SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad)258 IMPLICIT NONE259 260 c=======================================================================261 c262 c Calcul of equivalent solar angle and and fraction of day whithout263 c diurnal cycle.264 c265 c parmeters :266 c -----------267 c268 c Input :269 c -------270 c npts number of points271 c pdeclin solar declinaison272 c plat(npts) latitude273 c phaut hauteur typique de l'atmosphere274 c prad rayon planetaire '275 c276 c Output :277 c --------278 c pmu(npts) equivalent cosinus of the solar angle279 c pfract(npts) fractionnal day280 c281 c=======================================================================282 283 c-----------------------------------------------------------------------284 c285 c 0. Declarations :286 c -----------------287 288 c Arguments :289 c -----------290 INTEGER npts291 REAL plat(npts),pmu(npts), pfract(npts)292 REAL phaut,prad,pdeclin293 c294 c Local variables :295 c -----------------296 INTEGER j,i,ij,ig297 REAL pi298 REAL z,cz,sz,tz,phi,cphi,sphi,tphi299 REAL ap,a,t,b300 REAL alph301 302 c-----------------------------------------------------------------------303 304 !print*,'npts,pdeclin'305 !print*,npts,pdeclin306 pi = 4. * atan(1.0)307 !print*,'PI=',pi308 pi=2.*asin(1.)309 z = pdeclin310 cz = cos (z)311 sz = sin (z)312 !print*,'cz,sz',cz,sz313 314 ! DO j=jj_begin-offset,jj_end+offset315 ! DO i=ii_begin-offset,ii_end+offset316 ! ig=(j-1)*iim+i317 318 DO ig=1,npts319 phi = plat(ig)320 cphi = cos(phi)321 if (cphi.le.1.e-9) cphi=1.e-9322 sphi = sin(phi)323 tphi = sphi / cphi324 b = cphi * cz325 t = -tphi * sz / cz326 a = 1.0 - t*t327 ap = a328 IF(t.eq.0.) then329 t=0.5*pi330 ELSE331 IF (a.lt.0.) a = 0.332 t = sqrt(a) / t333 IF (t.lt.0.) then334 t = -atan (-t) + pi335 ELSE336 t = atan(t)337 ENDIF338 ENDIF339 340 pmu(ig) = (sphi*sz*t) / pi + b*sin(t)/pi341 pfract(ig) = t / pi342 IF (ap .lt.0.) then343 pmu(ig) = sphi * sz344 pfract(ig) = 1.0345 ENDIF346 IF (pmu(ig).le.0.0) pmu(ig) = 0.0347 pmu(ig) = pmu(ig) / pfract(ig)348 IF (pmu(ig).eq.0.) pfract(ig) = 0.349 ENDDO350 ! ENDDO351 c-----------------------------------------------------------------------352 c correction de rotondite:353 c ------------------------354 355 alph=phaut/prad356 ! DO j=jj_begin-offset,jj_end+offset357 ! DO i=ii_begin-offset,ii_end+offset358 ! ig=(j-1)*iim+i359 DO ig = 1,npts360 pmu(ig)=sqrt(1224.*pmu(ig)*pmu(ig)+1.)/35.361 362 c $ (sqrt(alph*alph*pmu(ig)*pmu(ig)+2.*alph+1.)-alph*pmu(ig))363 ENDDO364 ! ENDDO365 366 RETURN367 END SUBROUTINE mucorr368 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%369 370 SUBROUTINE sw(ngrid,nlayer,ldiurn,371 $ coefvis,albedo,372 $ plevel,ps_rad,pmu,pfract,psolarf0,373 $ fsrfvis,dtsw,374 $ lwrite)375 IMPLICIT NONE376 c=======================================================================377 c378 c Rayonnement solaire en atmosphere non diffusante avec un379 c coefficient d'absoprption gris.380 c'381 c=======================================================================382 c383 c declarations:384 c -------------385 c arguments:386 c ----------387 c388 INTEGER ngrid,nlayer,i,j,ij389 REAL albedo(ngrid),coefvis390 REAL pmu(ngrid),pfract(ngrid)391 REAL plevel(ngrid,nlayer+1),ps_rad392 REAL psolarf0393 REAL fsrfvis(ngrid),dtsw(ngrid,nlayer)394 LOGICAL lwrite,ldiurn395 c396 c variables locales:397 c ------------------398 c399 400 REAL zalb(ngrid),zmu(ngrid),zfract(ngrid)401 REAL zplev(ngrid,nlayer+1)402 REAL zflux(ngrid),zdtsw(ngrid,nlayer)403 404 INTEGER ig,l,nlevel,ncount,igout405 INTEGER,DIMENSION(ngrid)::index406 REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1)407 REAL zfsrfref(ngrid)408 REAL z1(ngrid)409 REAL zu(ngrid,nlayer+1)410 REAL tau0411 412 EXTERNAL SSUM413 EXTERNAL ismax414 REAL ismax415 416 LOGICAL firstcall417 SAVE firstcall418 DATA firstcall/.true./419 420 c-----------------------------------------------------------------------421 c 1. initialisations:422 c -------------------423 424 425 426 IF (firstcall) THEN427 428 IF (ngrid.NE.ngrid) THEN429 PRINT*,'STOP in inifis'430 PRINT*,'Probleme de dimenesions :'431 PRINT*,'ngrid = ',ngrid432 PRINT*,'ngrid = ',ngrid433 STOP434 ENDIF435 436 437 IF (nlayer.NE.llm) THEN438 PRINT*,'STOP in inifis'439 PRINT*,'Probleme de dimenesions :'440 PRINT*,'nlayer = ',nlayer441 PRINT*,'llm = ',llm442 STOP443 ENDIF444 445 ENDIF446 447 nlevel=nlayer+1448 c-----------------------------------------------------------------------449 c Definitions des tableaux locaux pour les points ensoleilles:450 c ------------------------------------------------------------451 452 IF (ldiurn) THEN453 ncount=0454 DO j=jj_begin-offset,jj_end+offset455 DO i=ii_begin-offset,ii_end+offset456 ig=(j-1)*iim+i457 index(ig)=0458 ENDDO459 ENDDO460 DO j=jj_begin-offset,jj_end+offset461 DO i=ii_begin-offset,ii_end+offset462 ig=(j-1)*iim+i463 IF(pfract(ig).GT.1.e-6) THEN464 ncount=ncount+1465 index(ncount)=ig466 ENDIF467 ENDDO468 ENDDO469 ! SARVESH CHANGED NCOUNT TO NGRID TO BE CONSISTENT ???470 471 CALL monGATHER(ncount,zfract,pfract,index)472 CALL monGATHER(ncount,zmu,pmu,index)473 CALL monGATHER(ncount,zalb,albedo,index)474 DO l=1,nlevel475 CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index)476 ENDDO477 ELSE478 ncount=ngrid479 zfract(:)=pfract(:)480 zmu(:)=pmu(:)481 zalb(:)=albedo(:)482 zplev(:,:)=plevel(:,:)483 ENDIF484 485 c-----------------------------------------------------------------------486 c calcul des profondeurs optiques integres depuis p=0:487 c ----------------------------------------------------488 489 tau0=-.5*log(coefvis)490 491 c calcul de la partie homogene de l'opacite492 c'493 tau0=tau0/ps_rad494 495 496 DO l=1,nlayer+1497 DO j=jj_begin-offset,jj_end+offset498 DO i=ii_begin-offset,ii_end+offset499 ig=(j-1)*iim+i500 zu(ig,l)=tau0*zplev(ig,l)501 ENDDO502 ENDDO503 ENDDO504 505 c-----------------------------------------------------------------------506 c 2. calcul de la transmission depuis le sommet de l'atmosphere:507 c' -----------------------------------------------------------508 509 DO l=1,nlevel510 DO j=jj_begin-offset,jj_end+offset511 DO i=ii_begin-offset,ii_end+offset512 ig=(j-1)*iim+i513 ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig))514 ENDDO515 ENDDO516 ENDDO517 518 IF (lwrite) THEN519 igout=ncount/2+1520 PRINT*521 PRINT*,'Diagnostique des transmission dans le spectre solaire'522 PRINT*,'zfract, zmu, zalb'523 PRINT*,zfract(igout), zmu(igout), zalb(igout)524 PRINT*,'Pression, quantite d abs, transmission'525 DO l=1,nlayer+1526 PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l)527 ENDDO528 ENDIF529 530 c-----------------------------------------------------------------------531 c 3. taux de chauffage, ray. solaire direct:532 c ------------------------------------------533 534 DO l=1,nlayer535 DO j=jj_begin-offset,jj_end+offset536 DO i=ii_begin-offset,ii_end+offset537 ig=(j-1)*iim+i538 zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)*539 $ (ztrdir(ig,l+1)-ztrdir(ig,l))/540 $ (cpp*(zplev(ig,l)-zplev(ig,l+1)))541 ENDDO542 ENDDO543 ENDDO544 IF (lwrite) THEN545 PRINT*546 PRINT*,'Diagnostique des taux de chauffage solaires:'547 PRINT*,' 1 taux de chauffage lie au ray. solaire direct'548 DO l=1,nlayer549 PRINT*,zdtsw(igout,l)550 ENDDO551 ENDIF552 553 554 c-----------------------------------------------------------------------555 c 4. calcul du flux solaire arrivant sur le sol:556 c ----------------------------------------------557 558 DO j=jj_begin-offset,jj_end+offset559 DO i=ii_begin-offset,ii_end+offset560 ig=(j-1)*iim+i561 z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1)562 zflux(ig)=(1.-zalb(ig))*z1(ig)563 zfsrfref(ig)= zalb(ig)*z1(ig)564 ENDDO565 ENDDO566 IF (lwrite) THEN567 PRINT*568 PRINT*,'Diagnostique des taux de chauffage solaires:'569 PRINT*,' 2 flux solaire net incident sur le sol'570 PRINT*,zflux(igout)571 ENDIF572 573 574 c-----------------------------------------------------------------------575 c 5.calcul des traansmissions depuis le sol, cas diffus:576 c ------------------------------------------------------577 578 DO l=1,nlevel579 DO j=jj_begin-offset,jj_end+offset580 DO i=ii_begin-offset,ii_end+offset581 ig=(j-1)*iim+i582 ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66)583 ENDDO584 ENDDO585 ENDDO586 587 IF (lwrite) THEN588 PRINT*589 PRINT*,'Diagnostique des taux de chauffage solaires'590 PRINT*,' 3 transmission avec les sol'591 PRINT*,'niveau transmission'592 DO l=1,nlevel593 PRINT*,l,ztrref(igout,l)594 ENDDO595 ENDIF596 597 c-----------------------------------------------------------------------598 c 6.ajout a l'echauffement de la contribution du ray. sol. reflechit:599 c' -------------------------------------------------------------------600 601 DO l=1,nlayer602 DO j=jj_begin-offset,jj_end+offset603 DO i=ii_begin-offset,ii_end+offset604 ig=(j-1)*iim+i605 zdtsw(ig,l)=zdtsw(ig,l)+606 $ g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/607 $ (cpp*(zplev(ig,l+1)-zplev(ig,l)))608 ENDDO609 ENDDO610 ENDDO611 612 c-----------------------------------------------------------------------613 c 10. sorites eventuelles:614 c ------------------------615 616 IF (lwrite) THEN617 PRINT*618 PRINT*,'Diagnostique des taux de chauffage solaires:'619 PRINT*,' 3 taux de chauffage total'620 DO l=1,nlayer621 PRINT*,zdtsw(igout,l)622 ENDDO623 ENDIF624 625 IF (ldiurn) THEN626 CALL zerophys(ngrid,fsrfvis)627 CALL monscatter(ncount,fsrfvis,index,zflux)628 CALL zerophys(ngrid*nlayer,dtsw)629 DO l=1,nlayer630 CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l))631 ENDDO632 ELSE633 !print*,'NOT DIURNE'634 fsrfvis(:)=zflux(:)635 dtsw(:,:)=zdtsw(:,:)636 ENDIF637 638 RETURN639 END SUBROUTINE sw640 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%641 SUBROUTINE lw(ngrid,nlayer,coefir,emissiv,642 $ pp,ps_rad,ptsurf,pt,643 $ pfluxir,pdtlw,644 $ lwrite)645 646 IMPLICIT NONE647 c=======================================================================648 c649 c calcul de l'evolution de la temperature sous l'effet du rayonnement650 c infra-rouge.651 c Pour simplifier, les transmissions sont precalculees et ne652 c dependent que de l'altitude.653 c654 c arguments:655 c ----------656 c'657 c entree:658 c -------659 c ngrid nombres de points de la grille horizontale660 c nlayer nombre de couches661 c ptsurf(ngrid) temperature de la surface662 c pt(ngrid,nlayer) temperature des couches663 c pp(ngrid,nlayer+1) pression entre les couches664 c lwrite variable logique pour sorties665 c666 c sortie:667 c -------668 c pdtlw(ngrid,nlayer) taux de refroidissement669 c pfluxir(ngrid) flux infrarouge sur le sol670 c671 c=======================================================================672 673 !c declarations:674 !c -------------675 !c arguments:676 !c' ----------677 678 INTEGER ngrid,nlayer679 REAL coefir,emissiv(ngrid),ps_rad680 REAL ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)681 REAL pdtlw(ngrid,nlayer),pfluxir(ngrid)682 LOGICAL lwrite683 684 c variables locales:685 c ------------------686 687 INTEGER nlevel,ilev,ig,i,il688 REAL zplanck(ngridmx,nlayermx+1),zcoef689 REAL zfluxup(ngridmx,nlayermx+1),zfluxdn(ngridmx,nlayermx+1)690 REAL zflux(ngridmx,nlayermx+1)691 REAL zlwtr1(ngridmx),zlwtr2(ngridmx)692 REAL zup(ngridmx,nlayermx+1),zdup(ngridmx)693 REAL stephan694 695 LOGICAL lstrong696 SAVE lstrong,stephan697 DATA lstrong/.true./698 c-----------------------------------------------------------------------699 c initialisations:700 c ----------------701 702 nlevel=nlayer+1703 stephan=5.67e-08704 705 706 pfluxir=0.0707 pdtlw=0.0708 !print*,"ngr,nlay,coefi",ngrid,nlayer,coefir709 c-----------------------------------------------------------------------710 c 2. calcul des quantites d'absorbants:711 c' -------------------------------------712 713 c absorption forte714 IF(lstrong) THEN715 DO ilev=1,nlevel716 DO ig=1,ngrid717 zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)718 ENDDO719 ENDDO720 IF(lwrite) THEN721 DO ilev=1,nlayer722 PRINT*,' up(',ilev,') = ',zup(ngrid/2+1,ilev)723 ENDDO724 ENDIF725 zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))726 727 c absorption faible728 ELSE729 DO ilev=1,nlevel730 DO ig=1,ngrid731 zup(ig,ilev)=pp(ig,ilev)732 ENDDO733 ENDDO734 zcoef=-log(coefir)/ps_rad735 ENDIF736 737 738 c-----------------------------------------------------------------------739 c 2. calcul de la fonction de corps noir:740 c ---------------------------------------741 742 DO ilev=1,nlayer743 DO ig=1,ngrid744 zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)745 zplanck(ig,ilev)=stephan*746 $ zplanck(ig,ilev)*zplanck(ig,ilev)747 ENDDO748 ENDDO749 750 c-----------------------------------------------------------------------751 c 4. flux descendants:752 c --------------------753 754 DO ilev=1,nlayer755 DO ig=1,ngrid756 zfluxdn(ig,ilev)=0.757 ENDDO758 DO ig=1,ngrid759 zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)760 ENDDO761 CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)762 763 DO il=nlayer,ilev,-1764 zlwtr2(:)=zlwtr1(:)765 DO ig=1,ngrid766 zdup(ig)=zup(ig,ilev)-zup(ig,il)767 ENDDO768 CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)769 DO ig=1,ngrid770 zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+771 $ zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))772 ENDDO773 ENDDO774 ENDDO775 776 DO ig=1,ngrid777 zfluxdn(ig,nlevel)=0.778 pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)779 ENDDO780 781 DO ig=1,ngrid782 zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)783 zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1)784 $ +(1.-emissiv(ig))*zfluxdn(ig,1)785 ENDDO786 787 c-----------------------------------------------------------------------788 c 3. flux montants:789 c ------------------790 791 DO ilev=1,nlayer792 DO ig=1,ngrid793 zdup(ig)=zup(ig,1)-zup(ig,ilev+1)794 ENDDO795 CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)796 DO ig=1,ngrid797 zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)798 ENDDO799 DO il=1,ilev800 zlwtr2(:)=zlwtr1(:)801 DO ig=1,ngrid802 zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)803 ENDDO804 CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)805 DO ig=1,ngrid806 zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+807 $ zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))808 ENDDO809 ENDDO810 811 ENDDO812 813 c-----------------------------------------------------------------------814 c 5. calcul des flux nets:815 c ------------------------816 817 DO ilev=1,nlevel818 DO ig=1,ngrid819 zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)820 ENDDO821 ENDDO822 823 c-----------------------------------------------------------------------824 c 6. Calcul des taux de refroidissement:825 c --------------------------------------826 827 DO ilev=1,nlayer828 DO ig=1,ngrid829 pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))*830 $ g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))831 ENDDO832 ENDDO833 834 c-----------------------------------------------------------------------835 c 10. sorties eventuelles:836 c ------------------------837 838 IF (lwrite) THEN839 840 PRINT*841 PRINT*,'Diagnostique rayonnement thermique'842 PRINT*843 PRINT*,'temperature ',844 $ 'flux montant flux desc. taux de refroid.'845 i=ngrid/2+1846 WRITE(6,9000) ptsurf(i)847 DO ilev=1,nlayer848 WRITE(6,'(i4,4e18.4)') ilev,pt(i,ilev),849 $ zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)850 ENDDO851 WRITE(6,9000) zfluxup(i,nlevel),zfluxdn(i,nlevel)852 853 ENDIF854 855 c-----------------------------------------------------------------------856 857 RETURN858 9000 FORMAT(4e18.4)859 END SUBROUTINE lw860 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%861 862 subroutine solang ( kgrid,psilon,pcolon,psilat,pcolat,863 & ptim1,ptim2,ptim3,pmu0,pfract )864 IMPLICIT NONE865 866 C867 C**** *LW* - ORGANIZES THE LONGWAVE CALCULATIONS868 C869 C PURPOSE.870 C --------871 C CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID872 C ==== INPUTS ===873 C874 C PSILON(KGRID) : SINUS OF THE LONGITUDE875 C PCOLON(KGRID) : COSINUS OF THE LONGITUDE876 C PSILAT(KGRID) : SINUS OF THE LATITUDE877 C PCOLAT(KGRID) : COSINUS OF THE LATITUDE878 C PTIM1 : SIN(DECLI)879 C PTIM2 : COS(DECLI)*COS(TIME)880 C PTIM3 : SIN(DECLI)*SIN(TIME)881 C882 C ==== OUTPUTS ===883 C884 C PMU0 (KGRID) : SOLAR ANGLE885 C PFRACT(KGRID) : DAY FRACTION OF THE TIME INTERVAL886 C887 C IMPLICIT ARGUMENTS : NONE888 C --------------------889 C890 C METHOD.891 C -------892 C893 C EXTERNALS.894 C ----------895 C896 C NONE897 C898 C REFERENCE.899 C ----------900 C901 C RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE902 C PALTRIDGE AND PLATT903 C904 C AUTHOR.905 C -------906 C FREDERIC HOURDIN907 C908 C MODIFICATIONS.909 C --------------910 C ORIGINAL :90-01-14911 C 92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher)912 C-----------------------------------------------------------------------913 C914 C ------------------------------------------------------------------915 916 C-----------------------------------------------------------------------917 C918 C* 0.1 ARGUMENTS919 C ---------920 C921 INTEGER kgrid922 REAL ptim1,ptim2,ptim3923 REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid)924 REAL psilat(kgrid), pcolat(kgrid)925 C926 INTEGER jl,i,j927 REAL ztim1,ztim2,ztim3928 C------------------------------------------------------------------------929 C------------------------------------------------------------------------930 C---------------------------------------------------------------------931 C932 C--------------------------------------------------------------------933 C934 C* 1. INITIALISATION935 C --------------936 C937 !!!!!! 100 CONTINUE938 C939 DO j=jj_begin-offset,jj_end+offset940 DO i=ii_begin-offset,ii_end+offset941 jl=(j-1)*iim+i942 pmu0(jl)=0.943 pfract(jl)=0.944 ENDDO945 ENDDO946 !C947 !C* 1.1 COMPUTATION OF THE SOLAR ANGLE948 !C ------------------------------949 !C950 DO j=jj_begin-offset,jj_end+offset951 DO i=ii_begin-offset,ii_end+offset952 jl=(j-1)*iim+i953 ztim1=psilat(jl)*ptim1954 ztim2=pcolat(jl)*ptim2955 ztim3=pcolat(jl)*ptim3956 pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl)957 ENDDO958 ENDDO959 !C960 !C* 1.2 DISTINCTION BETWEEN DAY AND NIGHT961 !C ---------------------------------962 !C963 DO j=jj_begin-offset,jj_end+offset964 DO i=ii_begin-offset,ii_end+offset965 jl=(j-1)*iim+i966 IF (pmu0(jl).gt.0.) THEN967 pfract(jl)=1.968 ELSE969 pmu0(jl)=0.970 pfract(jl)=0.971 ENDIF972 ENDDO973 ENDDO974 RETURN975 END SUBROUTINE solang976 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%977 978 SUBROUTINE monGATHER(n,a,b,index)979 c980 IMPLICIT NONE981 C982 INTEGER n,ng,index(n),i,j,ij983 REAL a(n),b(n)984 c985 DO 100 i=1,n986 a(i)=b(index(i))987 100 CONTINUE988 989 ! DO j=jj_begin-offset,jj_end+offset990 ! DO i=ii_begin-offset,ii_end+offset991 ! ij=(j-1)*iim+i992 ! a(ij)=b(index(ij))993 !c994 RETURN995 END SUBROUTINE monGATHER996 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%997 998 subroutine monscatter(n,a,index,b)999 C1000 implicit none1001 integer N,INDEX(n),I1002 real A(n),B(n)1003 c1004 c1005 DO 100 I=1,N1006 A(INDEX(I))=B(I)1007 100 CONTINUE1008 c1009 return1010 end SUBROUTINE monscatter1011 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1012 1013 SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)1014 IMPLICIT NONE1015 INTEGER ngrid1016 REAL coef1017 LOGICAL lstrong1018 REAL dup(ngrid),transm(ngrid)1019 INTEGER ig1020 1021 IF(lstrong) THEN1022 DO ig=1,ngrid1023 transm(ig)=exp(-coef*sqrt(dup(ig)))1024 ENDDO1025 ELSE1026 DO ig=1,ngrid1027 transm(ig)=exp(-coef*dup(ig))1028 ENDDO1029 ENDIF1030 RETURN1031 END subroutine lwtr1032 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1033 END MODULE RADIATION
Note: See TracChangeset
for help on using the changeset viewer.
