Ignore:
Timestamp:
Dec 20, 2019, 3:25:27 PM (5 years ago)
Author:
dubos
Message:

simple_physics : cleanup Mellor & Yamada

Location:
dynamico_lmdz/simple_physics
Files:
1 added
1 edited
1 moved

Legend:

Unmodified
Added
Removed

  • dynamico_lmdz/simple_physics/bash/replace_keywords.sh

    r4202 r4203  
    1010function replace_key()
    1111{
    12     LOW=$1
    13     UP=$2
     12    LOW="$1"
     13    UP="$2"
    1414    sed -i -e "s/\b${LOW}\b/${UP}/g" $3
    1515}
    … …  
    2222# enddo => ENDDO
    2323    replace_key subroutine SUBROUTINE $NEW
     24    replace_key 'do i=' 'DO i=' $NEW
     25    replace_key 'do k=' 'DO k=' $NEW
     26    replace_key 'end do' 'END DO' $NEW
     27    replace_key 'parameter' 'PARAMETER' $NEW
     28    replace_key 'integer' 'INTEGER' $NEW
     29    replace_key 'real' 'REAL' $NEW
     30    replace_key 'then' 'THEN' $NEW
     31    replace_key 'else' 'ELSE' $NEW
     32    replace_key 'endif' 'END IF' $NEW
     33    replace_key 'end if' 'END IF' $NEW
    2434    diff $NEW $ORIG
    2535    cp -i -u $NEW $ORIG

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

    r4202 r4203  
    1       subroutine my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh)
    2 
    3 ***********************************************************************
    4 ******* FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) *******
    5 ** q2 au interfaces entre mailles.
    6 ***********************************************************************
    7 
    8 
    9 ************** DECLARATIONS *******************************************
    10 
    11 
    12       integer imax,kmax
    13 
    14       parameter (impmax=5)
    15 
    16       real kappa,khmin,kmmin,kqmin,longc
    17 
    18       parameter (kappa=0.4)
    19       parameter (a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08,
    20      a           e1=1.8, e2=1.33)
    21       parameter (khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5,
    22      a           q2min=0.001, q2lmin=0.001)
    23       parameter (ghmax=0.023, ghmin=-0.28)
    24 
    25       real cd(imax),teta(imax,kmax),u(imax,kmax),v(imax,kmax),
    26      a     z(imax,kmax),zi(imax,kmax+1)
    27       real kh(imax,kmax+1),km(imax,kmax+1),q2(imax,kmax+1),
    28      a     long(imax,kmax+1)
    29       real unsdz(imax,kmax),unsdzi(imax,kmax+1)
    30       real kq(imax,kmax),
    31      a     m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1)
    32       real a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax),
    33      a     alph(imax,kmax)
    34       real ksdz2inf,ksdz2sup
    35 
    36 c      save q2,q2l
    37 
    38       save imp
    39 
    40       data imp/0/
    41 
    42 ***********************************************************************
    43 
    44       imp=imp+1
    45 
    46 ************** INCREMENTS VERTICAUX ***********************************
    47 
    48       do 9 i=1,imax
    49        zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax))
    50  9    continue
    51       do 10 k=1,kmax
    52        do 10 i=1,imax
    53         unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k))
    54  10   continue
    55       do 11 k=2,kmax
    56        do 11 i=1,imax
    57         unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1))
    58  11   continue
    59       do 12 i=1,imax
    60        unsdzi(i,1)=0.5/(z(i,1)-zi(i,1))
    61        unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax))
    62  12   continue
    63 
    64 ***********************************************************************
    65 
    66 
    67 ************** DIFFUSIVITES KH, KM et KQ ******************************
    68 * Ci-dessous, une premiere estimation des diffusivites turbulentes km *
    69 * et kh est effectuee pour utilisation dans les taux de production    *
    70 * et de destruction de q2 et q2l. On calcule aussi kq.                *
    71 
    72       do 100 k=2,kmax
    73        do 100 i=1,imax
    74         beta=2.0/(teta(i,k)+teta(i,k-1))
    75         n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1))
    76         n2(i,k)=amax1(0.0,n2(i,k))
    77         du=unsdzi(i,k)*(u(i,k)-u(i,k-1))
    78         dv=unsdzi(i,k)*(v(i,k)-v(i,k-1))
    79         m2(i,k)=du*du+dv*dv
    80         ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10)
    81         ri(i,k)=amax1(-0.1,min(4.0,ri(i,k)))
    82  100  continue
    83 
    84       do 110 k=2,kmax
    85        do 110 i=1,imax
    86         vt=kappa*(zi(i,k)-zi(i,1))
    87         long(i,k)=vt/(1.0+vt/160.0)
    88         if(n2(i,k).gt.0.0) then
    89          long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))
    90         endif
    91         gh=amax1(ghmin, 
    92      a     min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
    93         sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*
    94      a           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/
    95      a        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
    96         km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
    97         sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
    98         kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
    99  110  continue
    100       do 111 i=1,imax
    101        us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1)))
    102        vt1=(b1**0.666667)*us*us
    103        vt2=(b1**0.6666667)*kappa*kappa*
    104      a     m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1))
    105 c       q2(i,1)=amax1(q2min,vt1,vt2)
    106        q2(i,1)=amax1(q2min,vt1)
    107        long(i,1)=0.0
    108        long(i,kmax+1)=long(i,kmax)
    109        sq=0.2
    110        kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq)
    111  111  continue
    112 
    113       do 120 k=2,kmax
    114        do 120 i=1,imax
    115         longc=0.5*(long(i,k)+long(i,k+1))
    116         q2c=0.5*(q2(i,k)+q2(i,k+1))
    117         sq=0.2
    118         kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq)
    119  120  continue
    120 
    121 ***********************************************************************
    122 
    123 
    124 ************** CALCUL DE Q2 *******************************************
    125 
    126       do 200 k=2,kmax
    127        do 200 i=1,imax
    128         prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k)))
    129         dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k))
    130      a           +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k)))
    131         if(k.lt.kmax) then
    132          ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k)
    133         else
    134          ksdz2sup=0.0
    135         endif
    136         ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1)
    137         b(i,k)=-ksdz2inf*dt
    138         a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup)
    139         c(i,k)=-ksdz2sup*dt
    140         f(i,k)=q2(i,k)+dt*prod
    141  200  continue
    142       do 201 i=1,imax
    143        f(i,2)=f(i,2)
    144      a       +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1)
    145  201  continue
    146      
    147       do 210 i=1,imax
    148        alph(i,2)=a(i,2)
    149  210  continue
    150       do 211 k=3,kmax
    151        do 211 i=1,imax
    152         bet=b(i,k)/alph(i,k-1)
    153         alph(i,k)=a(i,k)-bet*c(i,k-1)
    154         f(i,k)=f(i,k)-bet*f(i,k-1)
    155  211  continue
    156      
    157       do 220 i=1,imax
    158        q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax))
    159        q2(i,kmax+1)=q2(i,kmax)
    160  220  continue
    161       do 221 k=kmax-1,2,-1
    162        do 221 i=1,imax
    163         q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k))
    164  221  continue
    165       do 222 i=1,imax
    166        q2(i,2)=amax1(q2(i,2),q2(i,1))
    167  222  continue
    168 
    169 ***********************************************************************
    170 
    171 
    172 ************** EVALUATION FINALE DE KH ET KM **************************
    173 
    174       do 400 k=2,kmax
    175        do 400 i=1,imax
    176         if(n2(i,k).gt.0.0) then
    177          long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))     
    178         endif
    179         gh=amax1(ghmin, 
    180      a         min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
    181         sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*
    182      a           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/
    183      a        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
    184         km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
    185         sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
    186         kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
    187  400  continue
    188       do 401 i=1,imax
    189        km(i,1)=kmmin
    190        km(i,kmax+1)=km(i,kmax)
    191        kh(i,1)=khmin
    192        kh(i,kmax+1)=kh(i,kmax)
    193  401  continue
    194  
    195 ***********************************************************************
    196 
    197 c      if(imp.eq.impmax) then
    198        am=1.0/float(imax)
    199        imp=0
    200        do 1000 k=kmax,1,-1
    201         au=0.0
    202         ateta=0.0
    203         az=0.0
    204         adz=0.0
    205         akq=0.0
    206         acd=0.0
    207         do 1001 i=1,imax
    208          au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k))
    209          ateta=ateta+am*teta(i,k)
    210          az=az+am*z(i,k)
    211          adz=adz+am*(1.0/unsdz(i,k))
    212          akq=akq+am*kq(i,k)
    213          acd=acd+am*cd(i)
    214  1001   continue
    215 c        write(*,2000) k,az,adz,au,ateta,akq,acd*1000.0
    216  2000   format(2x,i3,2x,6(f9.2,2x))
    217  1000  continue
    218        
    219        write(*,*)
    220        write(*,*)
    221 
    222        do 1002 k=kmax+1,1,-1
    223         azi=0.0
    224         adzi=0.0
    225         aq2=0.0
    226         al=0.0
    227         akm=0.0
    228         akh=0.0
    229         am2=0.0
    230         al0=0.0
    231         do 1003 i=1,imax
    232          azi=azi+am*zi(i,k)
    233          adzi=adzi+am*(1.0/unsdzi(i,k))
    234          aq2=aq2+am*q2(i,k)
    235          al=al+am*long(i,k)
    236          akm=akm+am*km(i,k)
    237          akh=akh+am*kh(i,k)
    238          am2=am2+am*m2(i,k)
    239 c         al0=al0+am*long0d(i)
    240  1003   continue
    241 c        write(*,2001) k,azi,aq2,al,akm,akh,am2*1.0e5
    242  2001   format(2x,i3,6(2x,f9.3))
    243  1002  continue
    244 c      endif
    245 
    246       return
    247       end
     1MODULE mellor_yamada
     2
     3  IMPLICIT NONE
     4
     5CONTAINS
     6
     7  SUBROUTINE my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh)
     8   
     9    !**********************************************************************
     10    !****** FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) *******
     11    !* q2 au interfaces entre mailles.                                     
     12    !**********************************************************************
     13   
     14   
     15    INTEGER, INTENT(IN) :: imax,kmax
     16    REAL, INTENT(IN) :: dt, gp
     17    REAL, INTENT(IN) :: z(imax,kmax), &
     18         u(imax,kmax), v(imax,kmax), teta(imax,kmax), cd(imax)         
     19    REAL, INTENT(INOUT) :: zi(imax,kmax+1), q2(imax,kmax+1)
     20    REAL, INTENT(OUT) :: long(imax,kmax+1), km(imax,kmax+1), kh(imax,kmax+1)
     21
     22    !************* DECLARATIONS *******************************************
     23   
     24    INTEGER, PARAMETER :: impmax=5
     25    REAL, PARAMETER :: kappa=0.4, &
     26         a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08,           &
     27         e1=1.8, e2=1.33, &                                       
     28         khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5,              &
     29         q2min=0.001, q2lmin=0.001, &                             
     30         ghmax=0.023, ghmin=-0.28
     31    REAL longc, au, ateta, az, adz, akq, acd, &
     32         adzi, aq2, al, akm, akh, am2, al0, &
     33         am, azi, bet, beta, dest, du, dv, gh, &
     34         prod, q2c, sh, sm, sq, us, vt, vt1, vt2
     35         
     36    REAL unsdz(imax,kmax),unsdzi(imax,kmax+1)
     37    REAL kq(imax,kmax),                                               &
     38         &     m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1)             
     39    REAL a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax),         &
     40         &     alph(imax,kmax)                                             
     41    REAL ksdz2inf,ksdz2sup
     42
     43    INTEGER :: i,k
     44   
     45    !************* INCREMENTS VERTICAUX ***********************************
     46   
     47    DO i=1,imax
     48       zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax))
     49    END DO
     50    DO k=1,kmax
     51       DO  i=1,imax
     52          unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k))
     53       END DO
     54    END DO
     55
     56    DO k=2,kmax
     57       DO i=1,imax
     58          unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1))
     59       END do
     60    END DO
     61
     62    DO i=1,imax
     63       unsdzi(i,1)=0.5/(z(i,1)-zi(i,1))
     64       unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax))
     65    END do
     66   
     67    !**********************************************************************
     68                     
     69   
     70    !************* DIFFUSIVITES KH, KM et KQ ******************************
     71    ! Ci-dessous, une premiere estimation des diffusivites turbulentes km *
     72    ! et kh est effectuee pour utilisation dans les taux de production    *
     73    ! et de destruction de q2 et q2l. On calcule aussi kq.                *
     74   
     75    DO k=2,kmax
     76       DO i=1,imax
     77          beta=2.0/(teta(i,k)+teta(i,k-1))
     78          n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1))
     79          n2(i,k)=amax1(0.0,n2(i,k))
     80          du=unsdzi(i,k)*(u(i,k)-u(i,k-1))
     81          dv=unsdzi(i,k)*(v(i,k)-v(i,k-1))
     82          m2(i,k)=du*du+dv*dv
     83          ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10)
     84          ri(i,k)=amax1(-0.1,min(4.0,ri(i,k)))
     85       END DO
     86    END DO
     87
     88    DO k=2,kmax
     89       DO i=1,imax
     90          vt=kappa*(zi(i,k)-zi(i,1))
     91          long(i,k)=vt/(1.0+vt/160.0)
     92          IF(n2(i,k).gt.0.0) THEN
     93             long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))
     94          END IF
     95          gh=amax1(ghmin,                                                 &
     96               &     min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))             
     97          sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*                          &
     98               &           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/     &
     99               &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))         
     100          km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
     101          sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
     102          kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
     103       END DO
     104    END DO
     105
     106    DO i=1,imax
     107       us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1)))
     108       vt1=(b1**0.666667)*us*us
     109       vt2=(b1**0.6666667)*kappa*kappa*                                 &
     110            &     m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1))                 
     111       !       q2(i,1)=amax1(q2min,vt1,vt2)                                   
     112       q2(i,1)=amax1(q2min,vt1)
     113       long(i,1)=0.0
     114       long(i,kmax+1)=long(i,kmax)
     115       sq=0.2
     116       kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq)
     117    END DO
     118   
     119    DO k=2,kmax
     120       DO i=1,imax
     121          longc=0.5*(long(i,k)+long(i,k+1))
     122          q2c=0.5*(q2(i,k)+q2(i,k+1))
     123          sq=0.2
     124          kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq)
     125       END DO
     126    END DO
     127
     128    !**********************************************************************
     129                                                                       
     130   
     131    !************* CALCUL DE Q2 *******************************************
     132   
     133    DO k=2,kmax
     134       DO i=1,imax
     135          prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k)))
     136          dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k))                            &
     137               &           +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k)))                 
     138          IF(k.lt.kmax) THEN
     139             ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k)
     140          ELSE
     141             ksdz2sup=0.0
     142          END IF
     143          ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1)
     144          b(i,k)=-ksdz2inf*dt
     145          a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup)
     146          c(i,k)=-ksdz2sup*dt
     147          f(i,k)=q2(i,k)+dt*prod
     148       END DO
     149    END DO
     150    DO i=1,imax
     151       f(i,2)=f(i,2)                                                    &
     152            &       +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1)                 
     153    END DO
     154   
     155    DO i=1,imax
     156       alph(i,2)=a(i,2)
     157    END DO
     158   
     159    DO k=3,kmax
     160       DO i=1,imax
     161          bet=b(i,k)/alph(i,k-1)
     162          alph(i,k)=a(i,k)-bet*c(i,k-1)
     163          f(i,k)=f(i,k)-bet*f(i,k-1)
     164       END DO
     165    END DO
     166   
     167    DO i=1,imax
     168       q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax))
     169       q2(i,kmax+1)=q2(i,kmax)
     170    END DO
     171   
     172    DO k=kmax-1,2,-1
     173       DO i=1,imax
     174          q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k))
     175       END DO
     176    END DO
     177   
     178    DO i=1,imax
     179       q2(i,2)=amax1(q2(i,2),q2(i,1))
     180    END DO
     181   
     182!**********************************************************************
     183                                                                       
     184                                                                       
     185!************* EVALUATION FINALE DE KH ET KM **************************
     186                                                                       
     187    DO k=2,kmax
     188       DO i=1,imax
     189          IF(n2(i,k).gt.0.0) THEN
     190             long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))   
     191          END IF
     192          gh=amax1(ghmin,                                                 &
     193               &         min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))         
     194          sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*                          &
     195               &           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/     &
     196               &        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))         
     197          km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
     198          sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
     199          kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
     200       END DO
     201    END DO
     202   
     203    DO i=1,imax
     204       km(i,1)=kmmin
     205       km(i,kmax+1)=km(i,kmax)
     206       kh(i,1)=khmin
     207       kh(i,kmax+1)=kh(i,kmax)
     208    END DO
     209   
     210    !**********************************************************************
     211   
     212    am=1.0/float(imax)
     213    DO k=kmax,1,-1
     214       au=0.0
     215       ateta=0.0
     216       az=0.0
     217       adz=0.0
     218       akq=0.0
     219       acd=0.0
     220       DO i=1,imax
     221          au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k))
     222          ateta=ateta+am*teta(i,k)
     223          az=az+am*z(i,k)
     224          adz=adz+am*(1.0/unsdz(i,k))
     225          akq=akq+am*kq(i,k)
     226          acd=acd+am*cd(i)
     227       END DO
     228    END DO
     229   
     230    DO k=kmax+1,1,-1
     231       azi=0.0
     232       adzi=0.0
     233       aq2=0.0
     234       al=0.0
     235       akm=0.0
     236       akh=0.0
     237       am2=0.0
     238       al0=0.0
     239       DO i=1,imax
     240          azi=azi+am*zi(i,k)
     241          adzi=adzi+am*(1.0/unsdzi(i,k))
     242          aq2=aq2+am*q2(i,k)
     243          al=al+am*long(i,k)
     244          akm=akm+am*km(i,k)
     245          akh=akh+am*kh(i,k)
     246          am2=am2+am*m2(i,k)
     247          !         al0=al0+am*long0d(i)                                         
     248       END DO
     249    END DO
     250   
     251  END SUBROUTINE my_25
     252 
     253END MODULE mellor_yamada
Note: See TracChangeset for help on using the changeset viewer.