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

simple_physics : cleanup convective adjustment

File:
1 moved

Legend:

Unmodified
Added
Removed

  • dynamico_lmdz/simple_physics/phyparam/param/convadj.F90

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