Changeset 4205


Ignore:
Timestamp:
Dec 20, 2019, 4:28:19 PM (5 years ago)
Author:
dubos
Message:

simple_physics : cleanup convective adjustment

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

Legend:

Unmodified
Added
Removed
  • dynamico_lmdz/simple_physics/phyparam/param/phyparam.F

    r4202 r4205  
    1414      USE astronomy
    1515      USE turbulence, ONLY : vdif
     16      USE convection, ONLY : convadj
    1617      USE solar, ONLY : solang, zenang, mucorr
    1718      USE radiative_sw, ONLY : sw
     
    166167
    167168
    168       EXTERNAL convadj
    169169      EXTERNAL ismin,ismax
    170170
  • dynamico_lmdz/simple_physics/phyparam/physics/convection.F90

    r4204 r4205  
    1 SUBROUTINE convadj(ngrid,nlay,ptimestep, &
    2      &                   pplay,pplev,ppopsk,   &
    3      &                  pu,pv,ph,              &
    4      &                  pdufi,pdvfi,pdhfi,     &
    5      &                  pduadj,pdvadj,pdhadj)
    6   USE phys_const
     1MODULE convection
    72  IMPLICIT NONE
    8  
    9   !=======================================================================
    10   !
    11   !   ajustement convectif sec
    12   !   on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement
    13   !
    14   !=======================================================================
    15  
    16   !-----------------------------------------------------------------------
    17   !   declarations:
    18   !   -------------
    19  
     3
     4CONTAINS
     5
     6  SUBROUTINE convadj(ngrid,nlay,ptimestep, &
     7       &                   pplay,pplev,ppopsk,   &
     8       &                  pu,pv,ph,              &
     9       &                  pdufi,pdvfi,pdhfi,     &
     10       &                  pduadj,pdvadj,pdhadj)
     11    USE phys_const
     12   
     13    !=======================================================================
     14    !
     15    !   ajustement convectif sec
     16    !   on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement
     17    !
     18    !=======================================================================
     19   
     20    !-----------------------------------------------------------------------
     21    !   declarations:
     22    !   -------------
     23   
    2024#include "dimensions.h"
    21  
    22   !   arguments:
    23   !   ----------
    24  
    25   INTEGER ngrid,nlay
    26   REAL ptimestep
    27   REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
    28   REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
    29   REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
    30   REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
    31  
    32   !   local:
    33   !   ------
    34  
    35   INTEGER ig,i,l,l1,l2,jj
    36   INTEGER jcnt, jadrs(ngrid)
    37  
    38   REAL*8 sig(nlay+1),sdsig(nlay),dsig(nlay)
    39   REAL*8 zu(ngrid,nlay),zv(ngrid,nlay)
    40   REAL*8 zh(ngrid,nlay)
    41   REAL*8 zu2(ngrid,nlay),zv2(ngrid,nlay)
    42   REAL*8 zh2(ngrid,nlay)
    43   REAL*8 zhm,zsm,zum,zvm,zalpha
    44  
    45   LOGICAL vtest(ngrid),down
    46  
    47   !
    48   !-----------------------------------------------------------------------
    49   !   initialisation:
    50   !   ---------------
    51   !
    52   !
    53   !-----------------------------------------------------------------------
    54   !   detection des profils a modifier:
    55   !   ---------------------------------
    56   !   si le profil est a modifier
    57   !   (i.e. ph(niv_sup) < ph(niv_inf) )
    58   !   alors le tableau "vtest" est mis a .TRUE. ;
    59   !   sinon, il reste a sa valeur initiale (.FALSE.)
    60   !   cette operation est vectorisable
    61   !   On en profite pour copier la valeur initiale de "ph"
    62   !   dans le champ de travail "zh"
    63  
    64   DO l=1,nlay
    65      DO ig=1,ngrid
    66         zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
    67         zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
    68         zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
    69      END DO
    70   END DO
    71  
    72   zu2(:,:)=zu(:,:)
    73   zv2(:,:)=zv(:,:)
    74   zh2(:,:)=zh(:,:)
    75  
    76   DO ig=1,ngrid
    77      vtest(ig)=.FALSE.
    78   END DO
    79   !
    80   DO l=2,nlay
    81      DO ig=1,ngrid
    82         !CRAY      vtest(ig)=CVMGM(.TRUE. , vtest(ig),
    83         !CRAY .                      zh2(ig,l)-zh2(ig,l-1))
    84         IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE.
    85      END DO
    86   END DO
    87   !
    88   !CRAY  CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)
    89   jcnt=0
    90   DO ig=1,ngrid
    91      IF(vtest(ig)) THEN
    92         jcnt=jcnt+1
    93         jadrs(jcnt)=ig
    94      ENDIF
    95   END DO
    96  
    97   !-----------------------------------------------------------------------
    98   !  Ajustement des "jcnt" profils instables indices par "jadrs":
    99   !  ------------------------------------------------------------
    100   !
    101   DO jj = 1, jcnt
    102      !
    103      i = jadrs(jj)
    104      !
    105      !   Calcul des niveaux sigma sur cette colonne
    106      DO l=1,nlay+1
    107         sig(l)=pplev(i,l)/pplev(i,1)
    108      ENDDO
    109      DO l=1,nlay
    110         dsig(l)=sig(l)-sig(l+1)
    111         sdsig(l)=ppopsk(i,l)*dsig(l)
    112      ENDDO
    113      l2 = 1
    114      !
    115      !      -- boucle de sondage vers le haut
    116      !
    117      DO WHILE(.TRUE.)
    118         ! 8000     CONTINUE
    119         !
    120         l2 = l2 + 1
    121         !
    122         IF (l2 .GT. nlay) EXIT ! Goto 8001
    123         !
    124         IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN
    125            !
    126            !          -- l2 est le niveau le plus haut de la colonne instable
    127            !
    128            l1 = l2 - 1
    129            l  = l1
    130            zsm = sdsig(l2)
    131            zhm = zh2(i, l2)
    132            !
    133            !          -- boucle de sondage vers le bas
    134            !
    135            DO WHILE(.TRUE.)                   
    136               ! 8020         CONTINUE
    137               zsm = zsm + sdsig(l)
    138               zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
    139               !
    140               !            -- doit on etendre la colonne vers le bas ?
    141               !
    142               !_EC (M1875) 20/6/87 : AND -> AND THEN
    143               !
    144               down = .FALSE.
    145               IF (l1 .NE. 1) THEN    !--  and then
    146                  IF (zhm .LT. zh2(i, l1-1)) THEN
    147                     down = .TRUE.
    148                  END IF
    149               END IF
    150               !
    151               IF (down) THEN
    152                  l1 = l1 - 1
    153                  l  = l1
    154               ELSE
    155                  !              -- peut on etendre la colonne vers le haut ?
    156                  IF (l2 .EQ. nlay) EXIT !Goto 8021
    157                  IF (zh2(i, l2+1) .GE. zhm) EXIT !Goto 8021
    158                  l2 = l2 + 1
    159                  l  = l2
    160               END IF
    161              
    162               ! GO TO 8020
    163            END DO
    164            ! 8021         CONTINUE
    165            !
    166            !          -- nouveau profil : constant (valeur moyenne)
    167            !
    168            zalpha=0.
    169            zum=0.
    170            zvm=0.
    171            DO l = l1, l2
    172               zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
    173               zh2(i, l) = zhm
    174               zum=zum+dsig(l)*zu(i,l)
    175               zvm=zvm+dsig(l)*zv(i,l)
    176            END DO
    177            zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
    178            zum=zum/(sig(l1)-sig(l2+1))
    179            zvm=zvm/(sig(l1)-sig(l2+1))
    180            IF(zalpha.GT.1.) THEN
    181               PRINT*,'WARNING dans convadj zalpha=',zalpha
    182               if(ig.eq.1) then
    183                  print*,'Au pole nord'
    184               elseif (ig.eq.ngrid) then
    185                  print*,'Au pole sud'
    186               else
    187                  print*,'Point i=', &
    188                       ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1
    189               endif
    190               STOP
    191               zalpha=1.
    192            ELSE
    193               !                IF(zalpha.LT.0.) STOP'zalpha=0'
    194               IF(zalpha.LT.1.e-5) zalpha=1.e-5
    195            ENDIF
    196            DO l=l1,l2
    197               zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
    198               zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
    199            ENDDO
    200            
    201            l2 = l2 + 1
    202            !
    203         END IF
    204         !
    205         !          GO TO 8000
    206      END DO
    207      ! 8001      CONTINUE
    208      !
    209      !
    210      DO l=1,nlay
    211         DO ig=1,ngrid
    212            pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
    213            pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
    214            pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
    215            !         pdhadj(ig,l)=0.
    216            !         pduadj(ig,l)=0.
    217            !         pdvadj(ig,l)=0.
    218         END DO
    219      END DO
    220      !
    221   END DO
     25   
     26    !   arguments:
     27    !   ----------
     28   
     29    INTEGER ngrid,nlay
     30    REAL ptimestep
     31    REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)
     32    REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)
     33    REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)
     34    REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)
     35   
     36    !   local:
     37    !   ------
     38   
     39    INTEGER ig,i,l,l1,l2,jj
     40    INTEGER jcnt, jadrs(ngrid)
     41   
     42    REAL*8 sig(nlay+1),sdsig(nlay),dsig(nlay)
     43    REAL*8 zu(ngrid,nlay),zv(ngrid,nlay)
     44    REAL*8 zh(ngrid,nlay)
     45    REAL*8 zu2(ngrid,nlay),zv2(ngrid,nlay)
     46    REAL*8 zh2(ngrid,nlay)
     47    REAL*8 zhm,zsm,zum,zvm,zalpha
     48   
     49    LOGICAL vtest(ngrid),down
     50   
     51    !
     52    !-----------------------------------------------------------------------
     53    !   initialisation:
     54    !   ---------------
     55    !
     56    !
     57    !-----------------------------------------------------------------------
     58    !   detection des profils a modifier:
     59    !   ---------------------------------
     60    !   si le profil est a modifier
     61    !   (i.e. ph(niv_sup) < ph(niv_inf) )
     62    !   alors le tableau "vtest" est mis a .TRUE. ;
     63    !   sinon, il reste a sa valeur initiale (.FALSE.)
     64    !   cette operation est vectorisable
     65    !   On en profite pour copier la valeur initiale de "ph"
     66    !   dans le champ de travail "zh"
     67   
     68    DO l=1,nlay
     69       DO ig=1,ngrid
     70          zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep
     71          zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep
     72          zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep
     73       END DO
     74    END DO
     75   
     76    zu2(:,:)=zu(:,:)
     77    zv2(:,:)=zv(:,:)
     78    zh2(:,:)=zh(:,:)
     79   
     80    DO ig=1,ngrid
     81       vtest(ig)=.FALSE.
     82    END DO
     83    !
     84    DO l=2,nlay
     85       DO ig=1,ngrid
     86          !CRAY      vtest(ig)=CVMGM(.TRUE. , vtest(ig),
     87          !CRAY .                      zh2(ig,l)-zh2(ig,l-1))
     88          IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE.
     89       END DO
     90    END DO
     91    !
     92    !CRAY  CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt)
     93    jcnt=0
     94    DO ig=1,ngrid
     95       IF(vtest(ig)) THEN
     96          jcnt=jcnt+1
     97          jadrs(jcnt)=ig
     98       ENDIF
     99    END DO
     100   
     101    !-----------------------------------------------------------------------
     102    !  Ajustement des "jcnt" profils instables indices par "jadrs":
     103    !  ------------------------------------------------------------
     104    !
     105    DO jj = 1, jcnt
     106       !
     107       i = jadrs(jj)
     108       !
     109       !   Calcul des niveaux sigma sur cette colonne
     110       DO l=1,nlay+1
     111          sig(l)=pplev(i,l)/pplev(i,1)
     112       ENDDO
     113       DO l=1,nlay
     114          dsig(l)=sig(l)-sig(l+1)
     115          sdsig(l)=ppopsk(i,l)*dsig(l)
     116       ENDDO
     117       l2 = 1
     118       !
     119       !      -- boucle de sondage vers le haut
     120       !
     121       DO WHILE(.TRUE.)
     122          ! 8000     CONTINUE
     123          !
     124          l2 = l2 + 1
     125          !
     126          IF (l2 .GT. nlay) EXIT ! Goto 8001
     127          !
     128          IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN
     129             !
     130             !          -- l2 est le niveau le plus haut de la colonne instable
     131             !
     132             l1 = l2 - 1
     133             l  = l1
     134             zsm = sdsig(l2)
     135             zhm = zh2(i, l2)
     136             !
     137             !          -- boucle de sondage vers le bas
     138             !
     139             DO WHILE(.TRUE.)                   
     140                ! 8020         CONTINUE
     141                zsm = zsm + sdsig(l)
     142                zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm
     143                !
     144                !            -- doit on etendre la colonne vers le bas ?
     145                !
     146                !_EC (M1875) 20/6/87 : AND -> AND THEN
     147                !
     148                down = .FALSE.
     149                IF (l1 .NE. 1) THEN    !--  and then
     150                   IF (zhm .LT. zh2(i, l1-1)) THEN
     151                      down = .TRUE.
     152                   END IF
     153                END IF
     154                !
     155                IF (down) THEN
     156                   l1 = l1 - 1
     157                   l  = l1
     158                ELSE
     159                   !              -- peut on etendre la colonne vers le haut ?
     160                   IF (l2 .EQ. nlay) EXIT !Goto 8021
     161                   IF (zh2(i, l2+1) .GE. zhm) EXIT !Goto 8021
     162                   l2 = l2 + 1
     163                   l  = l2
     164                END IF
     165               
     166                ! GO TO 8020
     167             END DO
     168             ! 8021         CONTINUE
     169             !
     170             !          -- nouveau profil : constant (valeur moyenne)
     171             !
     172             zalpha=0.
     173             zum=0.
     174             zvm=0.
     175             DO l = l1, l2
     176                zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l)
     177                zh2(i, l) = zhm
     178                zum=zum+dsig(l)*zu(i,l)
     179                zvm=zvm+dsig(l)*zv(i,l)
     180             END DO
     181             zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1)))
     182             zum=zum/(sig(l1)-sig(l2+1))
     183             zvm=zvm/(sig(l1)-sig(l2+1))
     184             IF(zalpha.GT.1.) THEN
     185                PRINT*,'WARNING dans convadj zalpha=',zalpha
     186                if(ig.eq.1) then
     187                   print*,'Au pole nord'
     188                elseif (ig.eq.ngrid) then
     189                   print*,'Au pole sud'
     190                else
     191                   print*,'Point i=', &
     192                        ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1
     193                endif
     194                STOP
     195                zalpha=1.
     196             ELSE
     197                !                IF(zalpha.LT.0.) STOP'zalpha=0'
     198                IF(zalpha.LT.1.e-5) zalpha=1.e-5
     199             ENDIF
     200             DO l=l1,l2
     201                zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l))
     202                zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l))
     203             ENDDO
     204             
     205             l2 = l2 + 1
     206             !
     207          END IF
     208          !
     209          !          GO TO 8000
     210       END DO
     211       ! 8001      CONTINUE
     212       !
     213       !
     214       DO l=1,nlay
     215          DO ig=1,ngrid
     216             pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep
     217             pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep
     218             pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep
     219             !         pdhadj(ig,l)=0.
     220             !         pduadj(ig,l)=0.
     221             !         pdvadj(ig,l)=0.
     222          END DO
     223       END DO
     224       !
     225    END DO
     226   
     227  END SUBROUTINE convadj
    222228
    223 END SUBROUTINE convadj
     229END MODULE convection
  • dynamico_lmdz/simple_physics/phyparam/physics/mellor_yamada.F90

    r4203 r4205  
    55CONTAINS
    66
    7   SUBROUTINE my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh)
     7  PURE SUBROUTINE my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh)
    88   
    99    !**********************************************************************
    1010    !****** FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) *******
    1111    !* q2 au interfaces entre mailles.                                     
    12     !**********************************************************************
    13    
     12    !**********************************************************************
    1413   
    1514    INTEGER, INTENT(IN) :: imax,kmax
     
    2928         q2min=0.001, q2lmin=0.001, &                             
    3029         ghmax=0.023, ghmin=-0.28
     30
    3131    REAL longc, au, ateta, az, adz, akq, acd, &
    3232         adzi, aq2, al, akm, akh, am2, al0, &
  • dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90

    r4201 r4205  
    1919    END DO
    2020  END SUBROUTINE monGATHER
     21
     22  PURE subroutine monscatter(n,a,index,b)
     23    INTEGER, INTENT(IN) :: n,index(n)
     24    REAL, INTENT(IN) :: b(n)
     25    REAL, INTENT(OUT) :: a(n)
     26    INTEGER :: i
     27    DO i=1,n
     28       a(index(i))=b(i)
     29    END DO
     30  end subroutine monscatter
    2131 
    2232  SUBROUTINE sw(ngrid,nlayer,ldiurn, &
     
    2535       fsrfvis,dtsw, &
    2636       lwrite)
    27     USE phys_const
     37    USE phys_const, ONLY : cpp, g
    2838    !=======================================================================
    2939    !
     
    242252 
    243253END MODULE radiative_sw
     254
Note: See TracChangeset for help on using the changeset viewer.