Changeset 4186


Ignore:
Timestamp:
Dec 19, 2019, 10:14:03 PM (5 years ago)
Author:
dubos
Message:

simple_physics : cleanup + bugfix vdif.F (wrong nb of args in CALL vdif_k)

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

Legend:

Unmodified
Added
Removed

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

    r4184 r4186  
    88      USE phys_const
    99      USE planet, ONLY : karman
     10      USE vdif_mod
    1011      IMPLICIT NONE
    1112
    … …  
    133134c
    134135      CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
    135 c     CALL my_25(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,ph,zcdv,
    136 c    a          pq2,pq2l,zkv,zkh)
    137         CALL vdif_k(ngrid,nlay,
    138      s   ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh,pq2,pq2l)
     136      CALL vdif_k(ngrid,nlay,
     137     s   ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh)
    139138
    140139      DO ig=1,ngrid

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

    r4184 r4186  
    1       SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
    2       USE planet, ONLY : karman
    3       IMPLICIT NONE
    4 c=======================================================================
    5 c
    6 c   Subject: computation of the surface drag coefficient using the
    7 c   -------  approch developed by Loui for ECMWF.
    8 c
    9 c   Author: Frederic Hourdin  15 /10 /93
    10 c   -------
    11 c
    12 c   Arguments:
    13 c   ----------
    14 c
    15 c   inputs:
    16 c   ------
    17 c     ngrid            size of the horizontal grid
    18 c     pg               gravity (m s -2)
    19 c     pz(ngrid)        height of the first atmospheric layer
    20 c     pu(ngrid)        u component of the wind in that layer
    21 c     pv(ngrid)        v component of the wind in that layer
    22 c     pts(ngrid)       surfacte temperature
    23 c     ph(ngrid)        potential temperature T*(p/ps)^kappa
    24 c
    25 c   outputs:
    26 c   --------
    27 c     pcdv(ngrid)      Cd for the wind
    28 c     pcdh(ngrid)      Cd for potential temperature
    29 c
    30 c=======================================================================
    31 c
    32 c-----------------------------------------------------------------------
    33 c   Declarations:
    34 c   -------------
     1MODULE vdif_mod
     2  IMPLICIT NONE
     3  SAVE
    354
    36 c   Arguments:
    37 c   ----------
     5CONTAINS
    386
    39       INTEGER ngrid,nlay
    40       REAL pz0(ngrid)
    41       REAL pg,pz(ngrid)
    42       REAL pu(ngrid),pv(ngrid)
    43       REAL pts(ngrid),ph(ngrid)
    44       REAL pcdv(ngrid),pcdh(ngrid)
    45 
    46 c   Local:
    47 c   ------
    48 
    49       INTEGER ig
    50 
    51       REAL zu2,z1,zri,zcd0,zz
    52 
    53       REAL b,c,d,c2b,c3bc,c3b,z0,umin2
    54       LOGICAL firstcal
    55       DATA b,c,d,umin2/5.,5.,5.,1.e-12/
    56       DATA firstcal/.true./
    57       SAVE b,c,d,c2b,c3bc,c3b,firstcal,umin2
    58 
    59 c-----------------------------------------------------------------------
    60 c   couche de surface:
    61 c   ------------------
    62 
    63 c     DO ig=1,ngrid
    64 c        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
    65 c        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
    66 c        pcdh(ig)=pcdv(ig)
    67 c     ENDDO
    68 c     RETURN
    69 
    70       IF (firstcal) THEN
    71          c2b=2.*b
    72          c3bc=3.*b*c
    73          c3b=3.*b
    74          firstcal=.false.
    75       ENDIF
    76 
    77 c!!!! WARNING, verifier la formule originale de Louis!
    78       DO ig=1,ngrid
    79          zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
    80          zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
    81          z1=1.+pz(ig)/pz0(ig)
    82          zcd0=karman/log(z1)
    83          zcd0=zcd0*zcd0*sqrt(zu2)
    84          IF(zri.LT.0.) THEN
    85             z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
    86             pcdv(ig)=zcd0*(1.-2.*z1)
    87             pcdh(ig)=zcd0*(1.-3.*z1)
    88          ELSE
    89             zz=sqrt(1.+d*zri)
    90             pcdv(ig)=zcd0/(1.+c2b*zri/zz)
    91             pcdh(ig)=zcd0/(1.+c3b*zri*zz)
    92          ENDIF
    93       ENDDO
    94 
    95 c-----------------------------------------------------------------------
    96 
    97       RETURN
    98       END
     7  SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
     8    USE planet, ONLY : karman
     9    !=======================================================================
     10    !
     11    !   Subject: computation of the surface drag coefficient using the
     12    !   -------  approch developed by Loui for ECMWF.
     13    !
     14    !   Author: Frederic Hourdin  15 /10 /93
     15    !   -------
     16    !
     17    !   Arguments:
     18    !   ----------
     19    !
     20    !   inputs:
     21    !   ------
     22    !     ngrid            size of the horizontal grid
     23    !     pg               gravity (m s -2)
     24    !     pz(ngrid)        height of the first atmospheric layer
     25    !     pu(ngrid)        u component of the wind in that layer
     26    !     pv(ngrid)        v component of the wind in that layer
     27    !     pts(ngrid)       surfacte temperature
     28    !     ph(ngrid)        potential temperature T*(p/ps)^kappa
     29    !
     30    !   outputs:
     31    !   --------
     32    !     pcdv(ngrid)      Cd for the wind
     33    !     pcdh(ngrid)      Cd for potential temperature
     34    !
     35    !=======================================================================
     36    !
     37    !-----------------------------------------------------------------------
     38    !   Declarations:
     39    !   -------------
     40   
     41    !   Arguments:
     42    !   ----------
     43   
     44    INTEGER ngrid,nlay
     45    REAL pz0(ngrid)
     46    REAL pg,pz(ngrid)
     47    REAL pu(ngrid),pv(ngrid)
     48    REAL pts(ngrid),ph(ngrid)
     49    REAL pcdv(ngrid),pcdh(ngrid)
     50   
     51    !   Local:
     52    !   ------
     53   
     54    INTEGER ig
     55   
     56    REAL zu2,z1,zri,zcd0,zz
     57   
     58    REAL b,c,d,c2b,c3bc,c3b,z0,umin2
     59    LOGICAL firstcal
     60    DATA b,c,d,umin2/5.,5.,5.,1.e-12/
     61    DATA firstcal/.true./
     62    SAVE b,c,d,c2b,c3bc,c3b,firstcal,umin2
     63   
     64    !-----------------------------------------------------------------------
     65    !   couche de surface:
     66    !   ------------------
     67   
     68    !     DO ig=1,ngrid
     69    !        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
     70    !        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
     71    !        pcdh(ig)=pcdv(ig)
     72    !     ENDDO
     73    !     RETURN
     74   
     75    IF (firstcal) THEN
     76       c2b=2.*b
     77       c3bc=3.*b*c
     78       c3b=3.*b
     79       firstcal=.false.
     80    ENDIF
     81   
     82!!!! WARNING, verifier la formule originale de Louis!
     83    DO ig=1,ngrid
     84       zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
     85       zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
     86       z1=1.+pz(ig)/pz0(ig)
     87       zcd0=karman/log(z1)
     88       zcd0=zcd0*zcd0*sqrt(zu2)
     89       IF(zri.LT.0.) THEN
     90          z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
     91          pcdv(ig)=zcd0*(1.-2.*z1)
     92          pcdh(ig)=zcd0*(1.-3.*z1)
     93       ELSE
     94          zz=sqrt(1.+d*zri)
     95          pcdv(ig)=zcd0/(1.+c2b*zri/zz)
     96          pcdh(ig)=zcd0/(1.+c3b*zri*zz)
     97       ENDIF
     98    ENDDO
     99   
     100    !-----------------------------------------------------------------------
     101   
     102    RETURN
     103  END SUBROUTINE vdif_cd
     104 
     105 
     106  SUBROUTINE vdif_k(ngrid,nlay,   &
     107       ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pcdv,pkv,pkh)
     108    ! FIXME : pkh :: pkv
     109    USE planet
     110    INTEGER ngrid,nlay
     111    REAL ptimestep
     112    REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
     113    REAL pz0(ngrid)
     114    REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
     115    REAL pg,pcdv(ngrid)
     116    REAL pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
     117   
     118    INTEGER ig,il
     119    REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix
     120   
     121    print*,'LMIXMIN',lmixmin
     122    DO ig=1,ngrid
     123       pkv(ig,1)=0.
     124       pkh(ig,1)=0.
     125       pkv(ig,nlay+1)=0.
     126       pkh(ig,nlay+1)=0.
     127    ENDDO
     128   
     129    DO il=2,nlay
     130       DO ig=1,ngrid
     131          z1=pzlev(ig,il)+pz0(ig)
     132          lmix=karman*z1/(1.+karman*z1/lmixmin)
     133          !           lmix=lmixmin
     134          ! WARNING test lmix=lmixmin
     135          zdu=pu(ig,il)-pu(ig,il-1)
     136          zdv=pv(ig,il)-pv(ig,il-1)
     137          zdz=pzlay(ig,il)-pzlay(ig,il-1)
     138          zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz)
     139          IF(zdvodz2.LT.1.e-5) THEN
     140             pkv(ig,il)=lmix*sqrt(emin_turb)
     141          ELSE
     142             zri=2.*pg*(ph(ig,il)-ph(ig,il-1)) &
     143                  / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2  )
     144             pkv(ig,il)= &
     145                  lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))
     146          ENDIF
     147          pkh(ig,il)=pkv(ig,il)
     148          !           IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
     149          !    s      zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
     150          !    s      lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
     151          !    s      ph(ig,il),ph(ig,il-1)
     152       ENDDO
     153    ENDDO
     154   
     155    RETURN
     156  END SUBROUTINE vdif_k
     157 
     158END MODULE vdif_mod
Note: See TracChangeset for help on using the changeset viewer.