Changeset 4199 for dynamico_lmdz


Ignore:
Timestamp:
Dec 20, 2019, 12:43:01 PM (5 years ago)
Author:
dubos
Message:

simple_physics : cleanup radiative transfer

Location:
dynamico_lmdz/simple_physics/phyparam
Files:
1 deleted
6 edited
1 moved

Legend:

Unmodified
Added
Removed

  • dynamico_lmdz/simple_physics/phyparam/param/phyparam.F

    r4198 r4199  
    1515      USE vdif_mod, ONLY : vdif
    1616      USE radiative, ONLY : mucorr, sw
     17      USE radiative_lw, ONLY : lw
     18
    1719c     
    1820      IMPLICIT NONE

  • dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90

    r4194 r4199  
    8585       WRITELOG(*,*) 'distance au soleil (en unite astronomique) :',pdist_sol
    8686       WRITELOG(*,*) 'declinaison (en degres) :',pdecli*180./pi
    87        CALL flush_log
     87       LOG_INFO('solarlong')
    8888    ENDIF
    8989   
    … …  
    168168    WRITELOG(*,*) 'longitude solaire du perihelie timeperi = ',timeperi
    169169   
    170     CALL flush_log
     170    LOG_INFO('iniorbit')
    171171   
    172172  END SUBROUTINE iniorbit

  • dynamico_lmdz/simple_physics/phyparam/physics/logging.F90

    r4194 r4199  
    1111
    1212  INTERFACE
    13      SUBROUTINE plugin(buf)
    14        CHARACTER(*), INTENT(IN) :: buf(:)
     13     SUBROUTINE plugin(lev, tag, buf)
     14       INTEGER, INTENT(IN) :: lev
     15       CHARACTER(*), INTENT(IN) :: tag, buf(:)
    1516     END SUBROUTINE plugin
    1617  END INTERFACE
    … …  
    2324  INTEGER :: logging_lineno=0
    2425
     26  INTEGER, PARAMETER, PUBLIC :: log_level_dbg=0, log_level_info=1
     27
    2528  PUBLIC :: logging_buf, logging_lineno, flush_log, flush_plugin
    2629
    2730CONTAINS
    2831
    29   SUBROUTINE flush_log
    30     CALL flush_plugin(logging_buf(1:logging_lineno))
     32  SUBROUTINE flush_log(lev,tag)
     33    INTEGER, INTENT(IN) :: lev
     34    CHARACTER(*), INTENT(IN) :: tag
     35    CALL flush_plugin(lev, TRIM(tag), logging_buf(1:logging_lineno))
    3136    logging_lineno=0
    3237  END SUBROUTINE flush_log
    3338
    34   SUBROUTINE default_flush_plugin(buf)
    35     CHARACTER(*), INTENT(IN) :: buf(:)
     39  SUBROUTINE default_flush_plugin(lev, tag, buf)
     40    INTEGER, INTENT(IN) :: lev
     41    CHARACTER(*), INTENT(IN) :: tag, buf(:)
    3642    INTEGER :: i
    3743    DO i=1, SIZE(buf)
    38        PRINT *, 'INFO : ', TRIM(buf(i))
     44       PRINT *, '[INFO ',tag,']', TRIM(buf(i))
    3945    END DO
    4046  END SUBROUTINE default_flush_plugin

  • dynamico_lmdz/simple_physics/phyparam/physics/radiative.F90

    r4198 r4199  
    55  PRIVATE
    66
    7   PUBLIC :: mucorr, sw
     7  PUBLIC :: mucorr, sw, lwtr
    88
    99CONTAINS
    … …  
    343343    END SUBROUTINE sw
    344344
     345    PURE SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)
     346      INTEGER, INTENT(IN) :: ngrid
     347      REAL,    INTENT(IN) :: coef
     348      LOGICAL, INTENT(IN) :: lstrong
     349      REAL,    INTENT(IN) :: dup(ngrid)
     350      REAL,    INTENT(OUT) :: transm(ngrid)
     351      INTEGER ig
     352      IF(lstrong) THEN
     353         DO ig=1,ngrid
     354            transm(ig)=exp(-coef*sqrt(dup(ig)))
     355         ENDDO
     356      ELSE
     357         DO ig=1,ngrid
     358            transm(ig)=exp(-coef*dup(ig))
     359         ENDDO
     360      ENDIF
     361     
     362    END SUBROUTINE lwtr
     363
    345364END MODULE radiative

  • dynamico_lmdz/simple_physics/phyparam/physics/radiative_lw.F90

    r4196 r4199  
    1       SUBROUTINE lw(ngrid,nlayer,coefir,emissiv,
    2      $             pp,ps_rad,ptsurf,pt,
    3      $             pfluxir,pdtlw,
    4      $             lwrite)
    5       USE phys_const
    6       IMPLICIT NONE
    7 c=======================================================================
    8 c
    9 c   calcul de l'evolution de la temperature sous l'effet du rayonnement
    10 c   infra-rouge.
    11 c   Pour simplifier, les transmissions sont precalculees et ne
    12 c   dependent que de l'altitude.
    13 c
    14 c   arguments:
    15 c   ----------
    16 c
    17 c   entree:
    18 c   -------
    19 c      ngrid             nombres de points de la grille horizontale
    20 c      nlayer              nombre de couches
    21 c      ptsurf(ngrid)     temperature de la surface
    22 c      pt(ngrid,nlayer)    temperature des couches
    23 c      pp(ngrid,nlayer+1)  pression entre les couches
    24 c      lwrite            variable logique pour sorties
    25 c
    26 c   sortie:
    27 c   -------
    28 c      pdtlw(ngrid,nlayer) taux de refroidissement
    29 c      pfluxir(ngrid)    flux infrarouge sur le sol
    30 c
    31 c=======================================================================
    32 
    33 c   declarations:
    34 c   -------------
    35 
    36 c   arguments:
    37 c   ----------
    38 
    39       INTEGER ngrid,nlayer
    40       REAL coefir,emissiv(ngrid),ps_rad
    41       REAL ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)
    42       REAL pdtlw(ngrid,nlayer),pfluxir(ngrid)
    43       LOGICAL lwrite
    44 
    45 c   variables locales:
    46 c   ------------------
    47 
    48       INTEGER nlevel,ilev,ig,i,il
    49       REAL zplanck(ngrid,nlayer+1),zcoef
    50       REAL zfluxup(ngrid,nlayer+1),zfluxdn(ngrid,nlayer+1)
    51       REAL zflux(ngrid,nlayer+1)
    52       REAL zlwtr1(ngrid),zlwtr2(ngrid)
    53       REAL zup(ngrid,nlayer+1),zdup(ngrid)
    54       REAL stephan
    55 
    56       LOGICAL lstrong
    57       SAVE lstrong,stephan
    58       DATA lstrong/.true./
    59 !$OMP THREADPRIVATE(lstrong,stephan)
    60 
    61 c-----------------------------------------------------------------------
    62 c   initialisations:
    63 c   ----------------
    64 
    65       nlevel=nlayer+1
    66       stephan=5.67e-08
    67 
    68 c-----------------------------------------------------------------------
    69 c   2. calcul des quantites d'absorbants:
    70 c   -------------------------------------
    71 
    72 c   absorption forte
    73       IF(lstrong) THEN
    74          DO ilev=1,nlevel
    75             DO ig=1,ngrid
    76                zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)
    77             ENDDO
    78          ENDDO
    79          IF(lwrite) THEN
    80             DO ilev=1,nlayer
    81              PRINT*,' up(',ilev,')  =  ',zup(ngrid/2+1,ilev)
    82             ENDDO
    83          ENDIF
    84          zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))
    85 
    86 c   absorption faible
    87       ELSE
    88          DO ilev=1,nlevel
    89             DO ig=1,ngrid
    90                zup(ig,ilev)=pp(ig,ilev)
    91             ENDDO
    92          ENDDO
    93          zcoef=-log(coefir)/ps_rad
    94       ENDIF
    95 
    96 
    97 c-----------------------------------------------------------------------
    98 c   2. calcul de la fonction de corps noir:
    99 c   ---------------------------------------
    100 
    101       DO ilev=1,nlayer
    102          DO ig=1,ngrid
    103             zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)
    104             zplanck(ig,ilev)=stephan*
    105      $      zplanck(ig,ilev)*zplanck(ig,ilev)
    106                  ENDDO
    107       ENDDO
    108 
    109 c-----------------------------------------------------------------------
    110 c   4. flux descendants:
    111 c   --------------------
    112 
    113       DO ilev=1,nlayer
    114          DO ig=1,ngrid
    115             zfluxdn(ig,ilev)=0.
    116                  ENDDO
    117          DO ig=1,ngrid
    118             zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)
    119          ENDDO
    120          CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
    121 
    122          DO il=nlayer,ilev,-1
    123             zlwtr2(:)=zlwtr1(:)
    124             DO ig=1,ngrid
    125                zdup(ig)=zup(ig,ilev)-zup(ig,il)
    126             ENDDO
    127             CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
    128             DO ig=1,ngrid
    129                zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+
    130      $         zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
    131                         ENDDO
    132                  ENDDO
    133       ENDDO
    134 
    135       DO ig=1,ngrid
    136          zfluxdn(ig,nlevel)=0.
    137          pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
    138       ENDDO
    139 
    140       DO ig=1,ngrid
    141          zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
    142          zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1)
    143      $   +(1.-emissiv(ig))*zfluxdn(ig,1)
    144       ENDDO
    145 
    146 c-----------------------------------------------------------------------
    147 c   3. flux montants:
    148 c   ------------------
    149 
    150       DO ilev=1,nlayer
    151          DO ig=1,ngrid
    152             zdup(ig)=zup(ig,1)-zup(ig,ilev+1)
    153          ENDDO
    154          CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
    155          DO ig=1,ngrid
    156             zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)
    157                  ENDDO
    158          DO il=1,ilev
    159             zlwtr2(:)=zlwtr1(:)
    160             DO ig=1,ngrid
    161                zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)
    162             ENDDO
    163             CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
    164             DO ig=1,ngrid
    165                zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+
    166      $         zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
    167                         ENDDO
    168          ENDDO
    169 
    170       ENDDO
    171 
    172 c-----------------------------------------------------------------------
    173 c   5. calcul des flux nets:
    174 c   ------------------------
    175 
    176       DO ilev=1,nlevel
    177          DO ig=1,ngrid
    178             zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)
    179                  ENDDO
    180       ENDDO
    181 
    182 c-----------------------------------------------------------------------
    183 c   6. Calcul des taux de refroidissement:
    184 c   --------------------------------------
    185    
    186       DO ilev=1,nlayer
    187          DO ig=1,ngrid
    188             pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))*
    189      $           g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))
    190                  ENDDO
    191       ENDDO
    192 
    193 c-----------------------------------------------------------------------
    194 c   10. sorties eventuelles:
    195 c   ------------------------
    196 
    197       IF (lwrite) THEN
    198 
    199          PRINT*
    200          PRINT*,'Diagnostique rayonnement thermique'
    201          PRINT*
    202          PRINT*,'temperature     ',
    203      $   'flux montant    flux desc.     taux de refroid.'
    204          i=ngrid/2+1
    205          WRITE(6,9000) ptsurf(i)
    206          DO ilev=1,nlayer
    207             WRITE(6,'(i4,4e18.4)') ilev,pt(i,ilev),
    208      $      zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
    209                  ENDDO
    210          WRITE(6,9000) zfluxup(i,nlevel),zfluxdn(i,nlevel)
    211 
    212       ENDIF
    213 
    214 c-----------------------------------------------------------------------
    215 
    216       RETURN
    217 9000  FORMAT(4e18.4)
    218       END
     1MODULE radiative_lw
     2
     3#include "use_logging.h"
     4
     5  IMPLICIT NONE
     6  SAVE
     7 
     8  PRIVATE
     9 
     10  PUBLIC :: lw
     11 
     12  LOGICAL, PARAMETER :: lstrong=.TRUE.
     13  REAL, PARAMETER :: stephan=5.67e-08
     14 
     15CONTAINS
     16 
     17  SUBROUTINE lw(ngrid,nlayer,coefir,emissiv, &
     18       pp,ps_rad,ptsurf,pt,              &
     19       pfluxir,pdtlw,                    &
     20       lwrite)
     21    USE phys_const, ONLY : cpp, g
     22    !=======================================================================
     23    !
     24    !   calcul de l'evolution de la temperature sous l'effet du rayonnement
     25    !   infra-rouge.
     26    !   Pour simplifier, les transmissions sont precalculees et ne
     27    !   dependent que de l'altitude.
     28    !
     29    !   arguments:
     30    !   ----------
     31    !
     32    !   entree:
     33    !   -------
     34    !      ngrid             nombres de points de la grille horizontale
     35    !      nlayer              nombre de couches
     36    !      ptsurf(ngrid)     temperature de la surface
     37    !      pt(ngrid,nlayer)    temperature des couches
     38    !      pp(ngrid,nlayer+1)  pression entre les couches
     39    !      lwrite            variable logique pour sorties
     40    !
     41    !   sortie:
     42    !   -------
     43    !      pdtlw(ngrid,nlayer) taux de refroidissement
     44    !      pfluxir(ngrid)    flux infrarouge sur le sol
     45    !
     46    !=======================================================================
     47   
     48    !   declarations:
     49    !   -------------
     50   
     51    !   arguments:
     52    !   ----------
     53   
     54    INTEGER, INTENT(IN)  ::  ngrid,nlayer
     55    REAL,    INTENT(IN)  :: coefir,emissiv(ngrid),ps_rad
     56    REAL,    INTENT(IN)  ::  ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1)
     57    REAL,    INTENT(OUT) ::  pdtlw(ngrid,nlayer),pfluxir(ngrid)
     58    LOGICAL, INTENT(IN)  :: lwrite
     59   
     60    !   variables locales:
     61    !   ------------------
     62   
     63    INTEGER nlevel,ilev,ig,i,il
     64    REAL zplanck(ngrid,nlayer+1),zcoef
     65    REAL zfluxup(ngrid,nlayer+1),zfluxdn(ngrid,nlayer+1)
     66    REAL zflux(ngrid,nlayer+1)
     67    REAL zlwtr1(ngrid),zlwtr2(ngrid)
     68    REAL zup(ngrid,nlayer+1),zdup(ngrid)
     69
     70    CHARACTER(:), PARAMETER :: tag='rad/lw'
     71    !-----------------------------------------------------------------------
     72    !   initialisations:
     73    !   ----------------
     74   
     75    nlevel=nlayer+1
     76   
     77    !-----------------------------------------------------------------------
     78    !   2. calcul des quantites d'absorbants:
     79    !   -------------------------------------
     80   
     81    !   absorption forte
     82    IF(lstrong) THEN
     83       DO ilev=1,nlevel
     84          DO ig=1,ngrid
     85             zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g)
     86          ENDDO
     87       ENDDO
     88       IF(lwrite) THEN
     89          DO ilev=1,nlayer
     90             WRITELOG(*,*) ' up(',ilev,')  =  ',zup(ngrid/2+1,ilev)
     91          ENDDO
     92          LOG_DBG(tag)
     93       ENDIF
     94       zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g))
     95       
     96       !   absorption faible
     97    ELSE
     98       DO ilev=1,nlevel
     99          DO ig=1,ngrid
     100             zup(ig,ilev)=pp(ig,ilev)
     101          ENDDO
     102       ENDDO
     103       zcoef=-log(coefir)/ps_rad
     104    ENDIF
     105   
     106   
     107    !-----------------------------------------------------------------------
     108    !   2. calcul de la fonction de corps noir:
     109    !   ---------------------------------------
     110   
     111    DO ilev=1,nlayer
     112       DO ig=1,ngrid
     113          zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev)
     114          zplanck(ig,ilev)=stephan* &
     115               zplanck(ig,ilev)*zplanck(ig,ilev)
     116       ENDDO
     117    ENDDO
     118   
     119    !-----------------------------------------------------------------------
     120    !   4. flux descendants:
     121    !   --------------------
     122   
     123    DO ilev=1,nlayer
     124       DO ig=1,ngrid
     125          zfluxdn(ig,ilev)=0.
     126       ENDDO
     127       DO ig=1,ngrid
     128          zdup(ig)=zup(ig,ilev)-zup(ig,nlevel)
     129       ENDDO
     130       CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
     131       
     132       DO il=nlayer,ilev,-1
     133          zlwtr2(:)=zlwtr1(:)
     134          DO ig=1,ngrid
     135             zdup(ig)=zup(ig,ilev)-zup(ig,il)
     136          ENDDO
     137          CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
     138          DO ig=1,ngrid
     139             zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+ &
     140                  zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
     141          ENDDO
     142       ENDDO
     143    ENDDO
     144   
     145    DO ig=1,ngrid
     146       zfluxdn(ig,nlevel)=0.
     147       pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1)
     148    ENDDO
     149   
     150    DO ig=1,ngrid
     151       zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig)
     152       zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1) &
     153            +(1.-emissiv(ig))*zfluxdn(ig,1)
     154    ENDDO
     155   
     156    !-----------------------------------------------------------------------
     157    !   3. flux montants:
     158    !   ------------------
     159   
     160    DO ilev=1,nlayer
     161       DO ig=1,ngrid
     162          zdup(ig)=zup(ig,1)-zup(ig,ilev+1)
     163       ENDDO
     164       CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
     165       DO ig=1,ngrid
     166          zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig)
     167       ENDDO
     168       DO il=1,ilev
     169          zlwtr2(:)=zlwtr1(:)
     170          DO ig=1,ngrid
     171             zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1)
     172          ENDDO
     173          CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1)
     174          DO ig=1,ngrid
     175             zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+ &
     176                  zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig))
     177          ENDDO
     178       ENDDO
     179       
     180    ENDDO
     181
     182    !-----------------------------------------------------------------------
     183    !   5. calcul des flux nets:
     184    !   ------------------------
     185   
     186    DO ilev=1,nlevel
     187       DO ig=1,ngrid
     188          zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev)
     189       ENDDO
     190    ENDDO
     191   
     192    !-----------------------------------------------------------------------
     193    !   6. Calcul des taux de refroidissement:
     194    !   --------------------------------------
     195   
     196    DO ilev=1,nlayer
     197       DO ig=1,ngrid
     198          pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))* &
     199               g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev)))
     200       ENDDO
     201    ENDDO
     202   
     203    !-----------------------------------------------------------------------
     204    !   10. sorties eventuelles:
     205    !   ------------------------
     206   
     207    IF (lwrite) THEN       
     208       WRITELOG(*,*) 'Diagnostique rayonnement thermique'
     209       WRITELOG(*,*) 'temperature     ', &
     210          'flux montant    flux desc.     taux de refroid.'
     211       i=ngrid/2+1
     212       WRITELOG(6,'(4e18.4)') ptsurf(i)
     213       DO ilev=1,nlayer
     214          WRITELOG(6,'(i4,4e18.4)') ilev,pt(i,ilev), &
     215               zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev)
     216       ENDDO
     217       WRITELOG(6,'(4e18.4)') zfluxup(i,nlevel),zfluxdn(i,nlevel)       
     218       LOG_DBG(tag)
     219    ENDIF
     220
     221!-----------------------------------------------------------------------
     222   
     223  END SUBROUTINE lw
     224
     225  PURE SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm)
     226    INTEGER, INTENT(IN) :: ngrid
     227    REAL,    INTENT(IN) :: coef
     228    LOGICAL, INTENT(IN) :: lstrong
     229    REAL,    INTENT(IN) :: dup(ngrid)
     230    REAL,    INTENT(OUT) :: transm(ngrid)
     231    INTEGER ig
     232    IF(lstrong) THEN
     233       DO ig=1,ngrid
     234          transm(ig)=exp(-coef*sqrt(dup(ig)))
     235       ENDDO
     236    ELSE
     237       DO ig=1,ngrid
     238          transm(ig)=exp(-coef*dup(ig))
     239       ENDDO
     240    ENDIF
     241   
     242  END SUBROUTINE lwtr
     243 
     244END MODULE radiative_lw

  • dynamico_lmdz/simple_physics/phyparam/physics/use_logging.h

    r4194 r4199  
    11   USE logging
    22#define WRITELOG(junk, fmt) logging_lineno = logging_lineno+1 ; WRITE(logging_buf(logging_lineno), fmt)
     3#define LOG_INFO(tag) CALL flush_log(log_level_info, tag)
     4#define LOG_DBG(tag) CALL flush_log(log_level_dbg, tag)

  • dynamico_lmdz/simple_physics/phyparam/physics/vdif_mod.F90

    r4194 r4199  
    262262               zb0(ig,ilev)
    263263       ENDDO
    264        CALL flush_log
     264       LOG_DBG('vdif')
    265265    ENDIF
    266266   
    … …  
    300300          PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
    301301       ENDDO
    302        CALL flush_log
     302       LOG_DBG('vdif')
    303303    ENDIF
    304304   
Note: See TracChangeset for help on using the changeset viewer.