Ignore:
Timestamp:
Mar 5, 2014, 2:19:12 PM (10 years ago)
Author:
lguez
Message:

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/yamada.F90

    r1988 r1992  
    1 !
     1
    22! $Header$
    3 !
    4       SUBROUTINE yamada(ngrid,dt,g,rconst,plev,temp
    5      s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,ustar
    6      s   ,l_mix)
    7       use dimphy
    8       IMPLICIT NONE
    9 c.......................................................................
    10 cym#include "dimensions.h"
    11 cym#include "dimphy.h"
    12 c.......................................................................
    13 c
    14 c dt : pas de temps
    15 c g  : g
    16 c zlev : altitude a chaque niveau (interface inferieure de la couche
    17 c        de meme indice)
    18 c zlay : altitude au centre de chaque couche
    19 c u,v : vitesse au centre de chaque couche
    20 c       (en entree : la valeur au debut du pas de temps)
    21 c teta : temperature potentielle au centre de chaque couche
    22 c        (en entree : la valeur au debut du pas de temps)
    23 c cd : cdrag
    24 c      (en entree : la valeur au debut du pas de temps)
    25 c q2 : $q^2$ au bas de chaque couche
    26 c      (en entree : la valeur au debut du pas de temps)
    27 c      (en sortie : la valeur a la fin du pas de temps)
    28 c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
    29 c      couche)
    30 c      (en sortie : la valeur a la fin du pas de temps)
    31 c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
    32 c      (en sortie : la valeur a la fin du pas de temps)
    33 c
    34 c.......................................................................
    35       REAL dt,g,rconst
    36       real plev(klon,klev+1),temp(klon,klev)
    37       real ustar(klon),snstable
    38       REAL zlev(klon,klev+1)
    39       REAL zlay(klon,klev)
    40       REAL u(klon,klev)
    41       REAL v(klon,klev)
    42       REAL teta(klon,klev)
    43       REAL cd(klon)
    44       REAL q2(klon,klev+1)
    45       REAL km(klon,klev+1)
    46       REAL kn(klon,klev+1)
    47       integer l_mix,ngrid
     3
     4SUBROUTINE yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
     5    cd, q2, km, kn, ustar, l_mix)
     6  USE dimphy
     7  IMPLICIT NONE
     8  ! .......................................................................
     9  ! ym#include "dimensions.h"
     10  ! ym#include "dimphy.h"
     11  ! .......................................................................
     12
     13  ! dt : pas de temps
     14  ! g  : g
     15  ! zlev : altitude a chaque niveau (interface inferieure de la couche
     16  ! de meme indice)
     17  ! zlay : altitude au centre de chaque couche
     18  ! u,v : vitesse au centre de chaque couche
     19  ! (en entree : la valeur au debut du pas de temps)
     20  ! teta : temperature potentielle au centre de chaque couche
     21  ! (en entree : la valeur au debut du pas de temps)
     22  ! cd : cdrag
     23  ! (en entree : la valeur au debut du pas de temps)
     24  ! q2 : $q^2$ au bas de chaque couche
     25  ! (en entree : la valeur au debut du pas de temps)
     26  ! (en sortie : la valeur a la fin du pas de temps)
     27  ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
     28  ! couche)
     29  ! (en sortie : la valeur a la fin du pas de temps)
     30  ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
     31  ! (en sortie : la valeur a la fin du pas de temps)
     32
     33  ! .......................................................................
     34  REAL dt, g, rconst
     35  REAL plev(klon, klev+1), temp(klon, klev)
     36  REAL ustar(klon), snstable
     37  REAL zlev(klon, klev+1)
     38  REAL zlay(klon, klev)
     39  REAL u(klon, klev)
     40  REAL v(klon, klev)
     41  REAL teta(klon, klev)
     42  REAL cd(klon)
     43  REAL q2(klon, klev+1)
     44  REAL km(klon, klev+1)
     45  REAL kn(klon, klev+1)
     46  INTEGER l_mix, ngrid
    4847
    4948
    50       integer nlay,nlev
    51 cym      PARAMETER (nlay=klev)
    52 cym      PARAMETER (nlev=klev+1)
     49  INTEGER nlay, nlev
     50  ! ym      PARAMETER (nlay=klev)
     51  ! ym      PARAMETER (nlev=klev+1)
    5352
    54       logical first
    55       save first
    56       data first/.true./
    57 c$OMP THREADPRIVATE(first)
     53  LOGICAL first
     54  SAVE first
     55  DATA first/.TRUE./
     56  !$OMP THREADPRIVATE(first)
    5857
    59       integer ig,k
     58  INTEGER ig, k
    6059
    61       real ri,zrif,zalpha,zsm
    62       real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
     60  REAL ri, zrif, zalpha, zsm
     61  REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)
    6362
    64       real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
    65       real l(klon,klev+1),l0(klon)
     63  REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
     64  REAL l(klon, klev+1), l0(klon)
    6665
    67       real sq(klon),sqz(klon),zz(klon,klev+1)
    68       integer iter
     66  REAL sq(klon), sqz(klon), zz(klon, klev+1)
     67  INTEGER iter
    6968
    70       real ric,rifc,b1,kap
    71       save ric,rifc,b1,kap
    72       data ric,rifc,b1,kap/0.195,0.191,16.6,0.3/
    73 c$OMP THREADPRIVATE(ric,rifc,b1,kap)
     69  REAL ric, rifc, b1, kap
     70  SAVE ric, rifc, b1, kap
     71  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.3/
     72  !$OMP THREADPRIVATE(ric,rifc,b1,kap)
    7473
    75       real frif,falpha,fsm
     74  REAL frif, falpha, fsm
    7675
    77       frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
    78       falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
    79       fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
     76  frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
     77  falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
     78  fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
    8079
    81       nlay=klev
    82       nlev=klev+1
    83      
    84       if (0.eq.1.and.first) then
    85       do ig=1,1000
    86          ri=(ig-800.)/500.
    87          if (ri.lt.ric) then
    88             zrif=frif(ri)
    89          else
    90             zrif=rifc
    91          endif
    92          if(zrif.lt.0.16) then
    93             zalpha=falpha(zrif)
    94             zsm=fsm(zrif)
    95          else
    96             zalpha=1.12
    97             zsm=0.085
    98          endif
    99          print*,ri,rif,zalpha,zsm
    100       enddo
    101       first=.false.
    102       endif
     80  nlay = klev
     81  nlev = klev + 1
    10382
    104 c  Correction d'un bug sauvage a verifier.
    105 c      do k=2,nlev
    106       do k=2,nlay
    107                                                           do ig=1,ngrid
    108          dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
    109          m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
    110      s             /(dz(ig,k)*dz(ig,k))
    111          n2(ig,k)=g*2.*(teta(ig,k)-teta(ig,k-1))
    112      s            /(teta(ig,k-1)+teta(ig,k))  /dz(ig,k)
    113          ri=n2(ig,k)/max(m2(ig,k),1.e-10)
    114          if (ri.lt.ric) then
    115             rif(ig,k)=frif(ri)
    116          else
    117             rif(ig,k)=rifc
    118          endif
    119          if(rif(ig,k).lt.0.16) then
    120             alpha(ig,k)=falpha(rif(ig,k))
    121             sm(ig,k)=fsm(rif(ig,k))
    122          else
    123             alpha(ig,k)=1.12
    124             sm(ig,k)=0.085
    125          endif
    126          zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
    127                                                           enddo
    128       enddo
     83  IF (0==1 .AND. first) THEN
     84    DO ig = 1, 1000
     85      ri = (ig-800.)/500.
     86      IF (ri<ric) THEN
     87        zrif = frif(ri)
     88      ELSE
     89        zrif = rifc
     90      END IF
     91      IF (zrif<0.16) THEN
     92        zalpha = falpha(zrif)
     93        zsm = fsm(zrif)
     94      ELSE
     95        zalpha = 1.12
     96        zsm = 0.085
     97      END IF
     98      PRINT *, ri, rif, zalpha, zsm
     99    END DO
     100    first = .FALSE.
     101  END IF
    129102
    130 c iterration pour determiner la longueur de melange
     103  ! Correction d'un bug sauvage a verifier.
     104  ! do k=2,nlev
     105  DO k = 2, nlay
     106    DO ig = 1, ngrid
     107      dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
     108      m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
     109        k-1))**2)/(dz(ig,k)*dz(ig,k))
     110      n2(ig, k) = g*2.*(teta(ig,k)-teta(ig,k-1))/(teta(ig,k-1)+teta(ig,k))/ &
     111        dz(ig, k)
     112      ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
     113      IF (ri<ric) THEN
     114        rif(ig, k) = frif(ri)
     115      ELSE
     116        rif(ig, k) = rifc
     117      END IF
     118      IF (rif(ig,k)<0.16) THEN
     119        alpha(ig, k) = falpha(rif(ig,k))
     120        sm(ig, k) = fsm(rif(ig,k))
     121      ELSE
     122        alpha(ig, k) = 1.12
     123        sm(ig, k) = 0.085
     124      END IF
     125      zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
     126    END DO
     127  END DO
    131128
    132                                                           do ig=1,ngrid
    133       l0(ig)=100.
    134                                                           enddo
    135       do k=2,klev-1
    136                                                           do ig=1,ngrid
    137         l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
    138                                                           enddo
    139       enddo
     129  ! iterration pour determiner la longueur de melange
    140130
    141       do iter=1,10
    142                                                           do ig=1,ngrid
    143          sq(ig)=1.e-10
    144          sqz(ig)=1.e-10
    145                                                           enddo
    146          do k=2,klev-1
    147                                                           do ig=1,ngrid
    148            q2(ig,k)=l(ig,k)**2*zz(ig,k)
    149            l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
    150      s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))
    151            zq=sqrt(q2(ig,k))
    152            sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
    153            sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
    154                                                           enddo
    155          enddo
    156                                                           do ig=1,ngrid
    157          l0(ig)=0.2*sqz(ig)/sq(ig)
    158                                                           enddo
    159 c(abd 3 5 2)         print*,'ITER=',iter,'  L0=',l0
     131  DO ig = 1, ngrid
     132    l0(ig) = 100.
     133  END DO
     134  DO k = 2, klev - 1
     135    DO ig = 1, ngrid
     136      l(ig, k) = l0(ig)*kap*zlev(ig, k)/(kap*zlev(ig,k)+l0(ig))
     137    END DO
     138  END DO
    160139
    161       enddo
     140  DO iter = 1, 10
     141    DO ig = 1, ngrid
     142      sq(ig) = 1.E-10
     143      sqz(ig) = 1.E-10
     144    END DO
     145    DO k = 2, klev - 1
     146      DO ig = 1, ngrid
     147        q2(ig, k) = l(ig, k)**2*zz(ig, k)
     148        l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
     149          k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
     150        zq = sqrt(q2(ig,k))
     151        sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
     152        sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
     153      END DO
     154    END DO
     155    DO ig = 1, ngrid
     156      l0(ig) = 0.2*sqz(ig)/sq(ig)
     157    END DO
     158    ! (abd 3 5 2)         print*,'ITER=',iter,'  L0=',l0
    162159
    163       do k=2,klev
    164                                                           do ig=1,ngrid
    165          l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
    166      s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))
    167          q2(ig,k)=l(ig,k)**2*zz(ig,k)
    168          km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
    169          kn(ig,k)=km(ig,k)*alpha(ig,k)
    170                                                           enddo
    171       enddo
     160  END DO
    172161
    173       return
    174       end
     162  DO k = 2, klev
     163    DO ig = 1, ngrid
     164      l(ig, k) = min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig, &
     165        k)+l0(ig)), 0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.E-10)))
     166      q2(ig, k) = l(ig, k)**2*zz(ig, k)
     167      km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
     168      kn(ig, k) = km(ig, k)*alpha(ig, k)
     169    END DO
     170  END DO
     171
     172  RETURN
     173END SUBROUTINE yamada
Note: See TracChangeset for help on using the changeset viewer.