Changeset 3757


Ignore:
Timestamp:
May 9, 2025, 9:22:41 AM (4 weeks ago)
Author:
emillour
Message:

Mars PCM:
More code tidying: puting routines in modules and modernizing some old
constructs. Tested to not change results with respect to previous version.
EM

Location:
trunk/LMDZ.MARS
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3756 r3757  
    48314831"blackl" routine (enforce "implicit none", make true constants "parameters")
    48324832and include it in the aerave module since it is only called there.
     4833
     4834== 09/05/2025 == EM
     4835More code tidying: puting routines in modules and modernizing some old
     4836constructs. Tested to not change results with respect to previous version.

  • trunk/LMDZ.MARS/libf/phymars/calldrag_noro_mod.F90

    r2642 r3757  
    3131      use dimradmars_mod, only: ndomainsz
    3232      use drag_noro_mod, only: drag_noro
     33      use sugwd_mod, only: sugwd
    3334      USE ioipsl_getin_p_mod, ONLY : getin_p
    3435     

  • trunk/LMDZ.MARS/libf/phymars/flusv.F

    r1268 r3757  
     1      MODULE flusv_mod
     2     
     3      IMPLICIT NONE
     4     
     5      CONTAINS
     6     
    17      SUBROUTINE flusv(KDLON,nsf,n,omega,g,tau,emis,bh,bsol,fah,fdh)
    2       use dimradmars_mod, only: ndlo2, ndlon, nflev
    3       IMPLICIT NONE
    4 c.......................................................................
    5 c
    6 c  calcul des flux ascendant et descendant aux interfaces entre n couches
    7 c  * dans l'infrarouge
    8 c  * B est une fonction lineaire de $\tau$ a l'interieur de chaque couche
    9 c  * le B du sol peut etre different de celui qui correspond au profil
    10 c    de la n-ieme couche
    11 c  * l'hypothese est une hypothese a deux flux isotropes sur chaque
    12 c    hemisphere ("hemispheric constant") + "source function technique"
    13 c    (voir Toon et al. 1988)
    14 c  * le flux descendant en haut de l'atmosphere est nul
    15 c  * les couches sont numerotees du haut de l'atmosphere vers le sol
    16 c
    17 c  in :   * KDLON      ---> dimension de vectorisation
     8        use dimradmars_mod, only: ndlo2, ndlon, nflev
     9        IMPLICIT NONE
     10c.......................................................................
     11c
     12c  compute the upward and downward fluxes at the interface between n layers
     13c  * in the infrared
     14c  * B is a linear function of $\tau$ in each layer
     15c  * B at the surface can be different than what corresponds to the profile
     16c    in the n-th layer
     17c  * the work hypothes isthat we have two isotropic fluxes for each
     18c    hemisphere ("hemispheric constant") + "technical source function"
     19c    (see Toon et al. 1988)
     20c  * the downward flux at the top of the atmosphere is zero
     21c  * layers are numbered from the top of the atmosphere to the ground
     22c
     23c  in :   * KDLON      ---> vectorisation dimension
    1824c         * nsf        ---> nsf=0 ==> "hemispheric constant"
    1925c                           nsf>0 ==> "hemispheric constant" + "source function"
    20 c         * n          ---> nombre de couches
    21 c         * omega(i)   ---> single scattering albedo pour la i-eme couche
    22 c         * g(i)       ---> asymmetry parameter pour la i-eme couche
    23 c         * tau(i)     ---> epaisseur optique de la i-eme couche
    24 c         * emis       ---> emissivite du sol
    25 c         * bh(i)      ---> luminance du corps noir en haut de la i-eme
    26 c                           couche, bh(n+1) pour la valeur au sol qui
    27 c                           correspond au profil de la n-ieme couche
    28 c         * bsol       ---> luminance du corps noir au sol
    29 c
    30 c  out :  * fah(i)     ---> flux ascendant en haut de la i-eme couche,
    31 c                           fah(n+1) pour le sol
    32 c         * fdh(i)     ---> flux descendant en haut de la i-eme couche,
    33 c                           fdh(n+1) pour le sol
    34 c
    35 c.......................................................................
    36 c  declaration des arguments
    37 c
    38       INTEGER KDLON,nsf,n
    39       REAL omega(NDLO2,n),g(NDLO2,n),tau(NDLO2,n),emis(NDLO2)
    40      &,bh(NDLO2,n+1),bsol(NDLO2),fah(NDLO2,n+1),fdh(NDLO2,n+1)
    41 c.......................................................................
    42 c  declaration des variables locales
    43 c
    44       REAL pi
    45       PARAMETER (pi=3.141592653589793E+0)
     26c         * n          ---> number of layers
     27c         * omega(i)   ---> single scattering albedo for the i-th layer
     28c         * g(i)       ---> asymmetry parameter for the i-th layer
     29c         * tau(i)     ---> optical thickness of the i-th layer
     30c         * emis       ---> ground emissivity
     31c         * bh(i)      ---> black body luminance at the top of the i-th layer,
     32c                           bh(n+1) for the ground value which
     33c                           corresponds to the profile for the n-th layer
     34c         * bsol       ---> black body luminance of the ground
     35c
     36c  out :  * fah(i)     ---> upward flux at the top of the i-th layer,
     37c                           fah(n+1) for the ground
     38c         * fdh(i)     ---> downward flux at the top of the i-th layer,
     39c                           fdh(n+1) for the ground
     40c
     41c.......................................................................
     42c  arguments
     43c
     44      INTEGER,INTENT(IN) :: KDLON,nsf,n
     45      REAL,INTENT(IN) :: omega(NDLO2,n),g(NDLO2,n)
     46      REAL,INTENT(IN) :: tau(NDLO2,n),emis(NDLO2)
     47      REAL,INTENT(IN) :: bh(NDLO2,n+1),bsol(NDLO2)
     48      REAL,INTENT(OUT) :: fah(NDLO2,n+1),fdh(NDLO2,n+1)
     49c.......................................................................
     50c  local variables
     51c
     52      REAL,PARAMETER :: pi=3.141592653589793E+0
    4653      INTEGER iv,i,j
    4754      REAL beta,gama1,gama2,amu1,grgama,b0,b1
    … …  
    5865     &    ,alpha1(NDLON,2*nflev),alpha2(NDLON,2*nflev)
    5966     &    ,sigma1(NDLON,2*nflev),sigma2(NDLON,2*nflev)
    60       INTEGER nq
    61       PARAMETER (nq=8)
    62       REAL x(nq),w(nq),gri(NDLON,nq)
    63       DATA x/1.9855071751231860E-2 , 0.1016667612931866E+0
     67      INTEGER,PARAMETER :: nq=8
     68      REAL,PARAMETER :: x(nq) =
     69     &     (/1.9855071751231860E-2 , 0.1016667612931866E+0
    6470     &     , 0.2372337950418355E+0 , 0.4082826787521751E+0
    6571     &     , 0.5917173212478250E+0 , 0.7627662049581645E+0
    66      &     , 0.8983332387068134E+0 , 0.9801449282487682E+0/
    67       DATA w/5.0614268145185310E-2 , 0.1111905172266872E+0
     72     &     , 0.8983332387068134E+0 , 0.9801449282487682E+0/)
     73      REAL,PARAMETER :: w(nq) =
     74     &     (/5.0614268145185310E-2 , 0.1111905172266872E+0
    6875     &     , 0.1568533229389437E+0 , 0.1813418916891810E+0
    6976     &     , 0.1813418916891810E+0 , 0.1568533229389437E+0
    70      &     , 0.1111905172266872E+0 , 5.0614268145185310E-2/
    71 c.......................................................................
    72 c
    73 c.......................................................................
    74       do 10001 i=1,n
    75                                                    do 99999 iv=1,KDLON
     77     &     , 0.1111905172266872E+0 , 5.0614268145185310E-2/)
     78      REAL :: gri(NDLON,nq)
     79c.......................................................................
     80c
     81c.......................................................................
     82      do i=1,n
     83        do iv=1,KDLON
    7684      beta=(1.E+0-g(iv,i))/2.E+0
    7785      gama1=2.E+0*(1.E+0-omega(iv,i)*(1.E+0-beta))
    … …  
    8189      grgama=(gama1-alambda(iv,i))/gama2
    8290c
    83 c ici une petite bidouille : si l'epaisseur optique d'une couche
    84 c est trop faible, $dB \over d\tau$ devient tres grand et le schema
    85 c s'ecroule. dans ces cas, on fait l'hypothese de couche isotherme.
     91c small hack here : if the optical depth of a layer is too small,
     92c $dB \over d\tau$ becomes very large and the scheme fails.
     93c In those cases we assume an isothermal layer.
    8694c
    8795      if (tau(iv,i).gt.1.E-3) then
    88       b0=bh(iv,i)
    89       b1=(bh(iv,i+1)-b0)/tau(iv,i)
     96        b0=bh(iv,i)
     97        b1=(bh(iv,i+1)-b0)/tau(iv,i)
    9098      else
    91       b0=(bh(iv,i)+bh(iv,i+1))/2.E+0
    92       b1=0.E+0
     99        b0=(bh(iv,i)+bh(iv,i+1))/2.E+0
     100        b1=0.E+0
    93101      endif
    94102c
    … …  
    111119      sigma2(iv,i)=alpha2(iv,i)
    112120c
    113 99999                                              continue
    114 10001 continue
    115 c.......................................................................
    116                                                    do 99998 iv=1,KDLON
    117       a(iv,1)=0.E+0
    118       b(iv,1)=e1(iv,1)
    119       d(iv,1)=-e2(iv,1)
    120       e(iv,1)=-cdh(iv,1)
    121 99998                                              continue
    122 c
    123       do 10002 i=1,n-1
    124       j=2*i+1
    125                                                    do 99997 iv=1,KDLON
    126       a(iv,j)=e2(iv,i)*e3(iv,i)-e4(iv,i)*e1(iv,i)
    127       b(iv,j)=e1(iv,i)*e1(iv,i+1)-e3(iv,i)*e3(iv,i+1)
    128       d(iv,j)=e3(iv,i)*e4(iv,i+1)-e1(iv,i)*e2(iv,i+1)
    129       e(iv,j)=e3(iv,i)*(cah(iv,i+1)-cab(iv,i))
    130      &       +e1(iv,i)*(cdb(iv,i)-cdh(iv,i+1))
    131 99997                                              continue
    132 10002 continue
    133 c
    134       do 10003 i=1,n-1
    135       j=2*i
    136                                                    do 99996 iv=1,KDLON
    137       a(iv,j)=e2(iv,i+1)*e1(iv,i)-e3(iv,i)*e4(iv,i+1)
    138       b(iv,j)=e2(iv,i)*e2(iv,i+1)-e4(iv,i)*e4(iv,i+1)
    139       d(iv,j)=e1(iv,i+1)*e4(iv,i+1)-e2(iv,i+1)*e3(iv,i+1)
    140       e(iv,j)=e2(iv,i+1)*(cah(iv,i+1)-cab(iv,i))
    141      &       +e4(iv,i+1)*(cdb(iv,i)-cdh(iv,i+1))
    142 99996                                              continue
    143 10003 continue
     121        enddo ! of do iv=1,KDLON
     122      enddo ! of do i=1,n
     123c.......................................................................
     124      do iv=1,KDLON
     125        a(iv,1)=0.E+0
     126        b(iv,1)=e1(iv,1)
     127        d(iv,1)=-e2(iv,1)
     128        e(iv,1)=-cdh(iv,1)
     129      enddo
     130c
     131      do i=1,n-1
     132        j=2*i+1
     133        do iv=1,KDLON
     134          a(iv,j)=e2(iv,i)*e3(iv,i)-e4(iv,i)*e1(iv,i)
     135          b(iv,j)=e1(iv,i)*e1(iv,i+1)-e3(iv,i)*e3(iv,i+1)
     136          d(iv,j)=e3(iv,i)*e4(iv,i+1)-e1(iv,i)*e2(iv,i+1)
     137          e(iv,j)=e3(iv,i)*(cah(iv,i+1)-cab(iv,i))
     138     &           +e1(iv,i)*(cdb(iv,i)-cdh(iv,i+1))
     139        enddo
     140      enddo ! of do i=1,n-1
     141c
     142      do i=1,n-1
     143        j=2*i
     144        do iv=1,KDLON
     145          a(iv,j)=e2(iv,i+1)*e1(iv,i)-e3(iv,i)*e4(iv,i+1)
     146          b(iv,j)=e2(iv,i)*e2(iv,i+1)-e4(iv,i)*e4(iv,i+1)
     147          d(iv,j)=e1(iv,i+1)*e4(iv,i+1)-e2(iv,i+1)*e3(iv,i+1)
     148          e(iv,j)=e2(iv,i+1)*(cah(iv,i+1)-cab(iv,i))
     149     &           +e4(iv,i+1)*(cdb(iv,i)-cdh(iv,i+1))
     150        enddo
     151      enddo ! of do i=1,n-1
    144152c     
    145153      j=2*n
    146                                                    do 99995 iv=1,KDLON
    147       a(iv,j)=e1(iv,n)-(1.E+0-emis(iv))*e3(iv,n)
    148       b(iv,j)=e2(iv,n)-(1.E+0-emis(iv))*e4(iv,n)
    149       d(iv,j)=0.E+0
    150       e(iv,j)=emis(iv)*pi*bsol(iv)-cab(iv,n)+(1.E+0-emis(iv))*cdb(iv,n)
    151 99995                                              continue
     154      do iv=1,KDLON
     155        a(iv,j)=e1(iv,n)-(1.E+0-emis(iv))*e3(iv,n)
     156        b(iv,j)=e2(iv,n)-(1.E+0-emis(iv))*e4(iv,n)
     157        d(iv,j)=0.E+0
     158        e(iv,j)=emis(iv)*pi*bsol(iv)-cab(iv,n)
     159     &          +(1.E+0-emis(iv))*cdb(iv,n)
     160      enddo
    152161c.......................................................................
    153162      call sys3v(KDLON,2*n,a,b,d,e,y)
    154163c.......................................................................
    155       do 10004 i=1,n
    156                                                    do 99994 iv=1,KDLON
    157       grg(iv,i)=grg(iv,i)*(y(iv,2*i-1)+y(iv,2*i))
    158       grh(iv,i)=grh(iv,i)*(y(iv,2*i-1)-y(iv,2*i))
    159       grj(iv,i)=grj(iv,i)*(y(iv,2*i-1)+y(iv,2*i))
    160       grk(iv,i)=grk(iv,i)*(y(iv,2*i-1)-y(iv,2*i))
    161 99994                                              continue
    162 10004 continue
    163 c.......................................................................
    164 c les valeurs de flux "hemispheric constant"
     164      do i=1,n
     165        do iv=1,KDLON
     166          grg(iv,i)=grg(iv,i)*(y(iv,2*i-1)+y(iv,2*i))
     167          grh(iv,i)=grh(iv,i)*(y(iv,2*i-1)-y(iv,2*i))
     168          grj(iv,i)=grj(iv,i)*(y(iv,2*i-1)+y(iv,2*i))
     169          grk(iv,i)=grk(iv,i)*(y(iv,2*i-1)-y(iv,2*i))
     170        enddo
     171      enddo
     172c.......................................................................
     173c values of "hemispheric constant" fluxes
    165174c
    166175      IF (nsf.eq.0) THEN
    167       do 10005 i=1,n
    168                                                    do 99993 iv=1,KDLON
    169       fah(iv,i)=e3(iv,i)*y(iv,2*i-1)-e4(iv,i)*y(iv,2*i)+cah(iv,i)
    170       fdh(iv,i)=e1(iv,i)*y(iv,2*i-1)-e2(iv,i)*y(iv,2*i)+cdh(iv,i)
    171 99993                                              continue
    172 10005 continue
    173                                                    do 99992 iv=1,KDLON
    174       fah(iv,n+1)=e1(iv,n)*y(iv,2*n-1)+e2(iv,n)*y(iv,2*n)+cab(iv,n)
    175       fdh(iv,n+1)=e3(iv,n)*y(iv,2*n-1)+e4(iv,n)*y(iv,2*n)+cdb(iv,n)
    176 99992                                              continue
    177       GOTO 10100
    178       ENDIF
    179 c.......................................................................
    180 c passage a "source function"
    181 c
    182 c on applique une quadrature de dimension nq
    183 c x est le vecteur des \mu de la quadrature
    184 c w est le vecteur des poids correspondants
    185 c nq est fixe en parameter
    186 c x et w sont fixes en data
    187 c
    188 c.......................................................................
    189 c on part d'en haut et on descent selon les nq angles pour calculer
    190 c tous les flux descendants
    191 c
    192       do 10006 j=1,nq
    193                                                    do 99991 iv=1,KDLON
    194       gri(iv,j)=0.E+0
    195 99991                                              continue
    196 10006 continue
    197                                                    do 99990 iv=1,KDLON
    198       fdh(iv,1)=0.E+0
    199 99990                                              continue
    200       do 10007 i=1,n
    201       do 10008 j=1,nq
    202                                                    do 99989 iv=1,KDLON
    203       gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j))
    204      &+grj(iv,i)/(alambda(iv,i)*x(j)+1.E+0)
    205      &*(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j))))
    206      &+grk(iv,i)/(alambda(iv,i)*x(j)-1.E+0)
    207      &*(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i)))
    208      &+sigma1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j)))
    209      &+sigma2(iv,i)*(x(j)*exp(-tau(iv,i)/x(j))+tau(iv,i)-x(j))
    210 99989                                              continue
    211 10008 continue
    212                                                    do 99988 iv=1,KDLON
    213       fdh(iv,i+1)=0.E+0
    214 99988                                              continue
    215       do 10009 j=1,nq
    216                                                    do 99987 iv=1,KDLON
    217       fdh(iv,i+1)=fdh(iv,i+1)+w(j)*x(j)*gri(iv,j)
    218 99987                                              continue
    219 10009 continue
    220 10007 continue
    221 c.......................................................................
    222 c on applique la condition de reflexion a sol
    223 c
    224                                                    do 99986 iv=1,KDLON
    225       fah(iv,n+1)=(1.E+0-emis(iv))*fdh(iv,n+1)+pi*emis(iv)*bsol(iv)
    226 99986                                              continue
    227       do 10010 j=1,nq
    228                                                    do 99985 iv=1,KDLON
    229       gri(iv,j)=2.E+0*fah(iv,n+1)
    230 99985                                              continue
    231 10010 continue
    232 c.......................................................................
    233 c on remonte pour calculer tous les flux ascendants
     176        do i=1,n
     177          do iv=1,KDLON
     178            fah(iv,i)=e3(iv,i)*y(iv,2*i-1)-e4(iv,i)*y(iv,2*i)+cah(iv,i)
     179            fdh(iv,i)=e1(iv,i)*y(iv,2*i-1)-e2(iv,i)*y(iv,2*i)+cdh(iv,i)
     180          enddo
     181        enddo
     182        do iv=1,KDLON
     183          fah(iv,n+1)=e1(iv,n)*y(iv,2*n-1)+e2(iv,n)*y(iv,2*n)+cab(iv,n)
     184          fdh(iv,n+1)=e3(iv,n)*y(iv,2*n-1)+e4(iv,n)*y(iv,2*n)+cdb(iv,n)
     185        enddo
     186      ELSE
     187c.......................................................................
     188c going to the "source function"
     189c
     190c apply a quadrature over nq (fixed parameter) points
     191c x is the vector of the \mu of the quadrature
     192c w is the vector of corresponding weights
     193c x() et w() are fixed parameters
     194c
     195c.......................................................................
     196c start from the top and go down along the nq angles to compute all
     197c downward fluxes
     198c
     199      do j=1,nq
     200        do iv=1,KDLON
     201          gri(iv,j)=0.E+0
     202        enddo
     203      enddo
     204      do iv=1,KDLON
     205        fdh(iv,1)=0.E+0
     206      enddo
     207      do i=1,n
     208        do j=1,nq
     209          do iv=1,KDLON
     210            gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j))
     211     &      +grj(iv,i)/(alambda(iv,i)*x(j)+1.E+0)
     212     &      *(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j))))
     213     &      +grk(iv,i)/(alambda(iv,i)*x(j)-1.E+0)
     214     &      *(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i)))
     215     &      +sigma1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j)))
     216     &      +sigma2(iv,i)*(x(j)*exp(-tau(iv,i)/x(j))+tau(iv,i)-x(j))
     217          enddo ! of do iv=1,KDLON
     218        enddo ! of do j=1,nq
     219        do iv=1,KDLON
     220          fdh(iv,i+1)=0.E+0
     221        enddo
     222        do j=1,nq
     223          do iv=1,KDLON
     224            fdh(iv,i+1)=fdh(iv,i+1)+w(j)*x(j)*gri(iv,j)
     225          enddo ! of do iv=1,KDLON
     226        enddo ! of do j=1,nq
     227      enddo ! of do i=1,n
     228c.......................................................................
     229c apply the reflexion condition on the ground
     230c
     231      do iv=1,KDLON
     232        fah(iv,n+1)=(1.E+0-emis(iv))*fdh(iv,n+1)+pi*emis(iv)*bsol(iv)
     233      enddo
     234      do j=1,nq
     235        do iv=1,KDLON
     236          gri(iv,j)=2.E+0*fah(iv,n+1)
     237        enddo
     238      enddo
     239c.......................................................................
     240c going back up to compute all the upward fluxes
    234241c 
    235       do 10011 i=n,1,-1
    236       do 10012 j=1,nq
    237                                                    do 99984 iv=1,KDLON
    238       gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j))
    239      &+grg(iv,i)/(alambda(iv,i)*x(j)-1.E+0)
    240      &*(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i)))
    241      &+grh(iv,i)/(alambda(iv,i)*x(j)+1.E+0)
    242      &*(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j))))
    243      &+alpha1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j)))
    244      &+alpha2(iv,i)*(x(j)-(tau(iv,i)+x(j))*exp(-tau(iv,i)/x(j)))
    245 99984                                              continue
    246 10012 continue
    247                                                    do 99983 iv=1,KDLON
    248       fah(iv,i)=0.E+0
    249 99983                                              continue
    250       do 10013 j=1,nq
    251                                                    do 99982 iv=1,KDLON
    252       fah(iv,i)=fah(iv,i)+w(j)*x(j)*gri(iv,j)
    253 99982                                              continue
    254 10013 continue
    255 10011 continue
    256 c.......................................................................
    257 10100 continue     
    258 c.......................................................................
    259 c
    260 c.......................................................................
    261       RETURN
    262       END
     242      do i=n,1,-1
     243        do j=1,nq
     244          do iv=1,KDLON
     245            gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j))
     246     &      +grg(iv,i)/(alambda(iv,i)*x(j)-1.E+0)
     247     &      *(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i)))
     248     &      +grh(iv,i)/(alambda(iv,i)*x(j)+1.E+0)
     249     &      *(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j))))
     250     &      +alpha1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j)))
     251     &      +alpha2(iv,i)*(x(j)-(tau(iv,i)+x(j))*exp(-tau(iv,i)/x(j)))
     252          enddo ! of do iv=1,KDLON
     253        enddo ! of do j=1,nq
     254        do iv=1,KDLON
     255          fah(iv,i)=0.E+0
     256        enddo
     257        do j=1,nq
     258          do iv=1,KDLON
     259            fah(iv,i)=fah(iv,i)+w(j)*x(j)*gri(iv,j)
     260          enddo
     261        enddo ! of do j=1,nq
     262      enddo ! of do i=n,1,-1
     263c.......................................................................
     264      ENDIF ! of IF (nsf.eq.0)
     265c.......................................................................
     266c
     267c.......................................................................
     268
     269      END SUBROUTINE flusv
    263270
    264271c ***************************************************************
    … …  
    270277c.......................................................................
    271278c
    272 c  inversion d'un systeme lineaire tridiagonal
     279c  solve a tridiagonal linear system such that:
    273280c
    274281c  |  b1  d1               |   | y1 |   | e1 |
    … …  
    278285c  |              an  bn   |   | yn |   | en |
    279286c
    280 c  in :   * KDLON        --> dimension de vectorisation
    281 c         * n            --> dimension du systeme
    282 c         * a,b,d,e      --> voir dessin
    283 c
    284 c  out :  * y            --> voir dessin
    285 c
    286 c.......................................................................
    287 c  declaration des arguments
    288 c
    289       INTEGER KDLON,n
    290       REAL a(NDLO2,n),b(NDLO2,n),d(NDLO2,n),e(NDLO2,n),y(NDLO2,n)
    291 c.......................................................................
    292 c  declaration des variables locales
    293 c
    294       INTEGER iv,i
    295       REAL as(NDLON,4*nflev),ds(NDLON,4*nflev)
     287c  in :   * KDLON        --> vectorisation dimension
     288c         * n            --> system size
     289c         * a,b,d,e      --> coefficients as shown above
     290c
     291c  out :  * y            --> see above
     292c
     293c.......................................................................
     294c  arguments
     295c
     296      INTEGER,INTENT(IN) :: KDLON,n
     297      REAL,INTENT(IN) :: a(NDLO2,n),b(NDLO2,n),d(NDLO2,n),e(NDLO2,n)
     298      REAL,INTENT(OUT) :: y(NDLO2,n)
     299c.......................................................................
     300c  local variables
     301c
     302      INTEGER :: iv,i
     303      REAL :: as(NDLON,4*nflev),ds(NDLON,4*nflev)
    296304     &    ,x(NDLON,4*nflev)
    297305c.......................................................................
    298306c
    299307c.......................................................................
    300                                                    do 99999 iv=1,KDLON
    301       as(iv,n)=a(iv,n)/b(iv,n)
    302       ds(iv,n)=e(iv,n)/b(iv,n)
    303 99999                                              continue
    304       do 10001 i=n-1,1,-1
    305                                                    do 99998 iv=1,KDLON
    306       x(iv,i)=1.E+0/(b(iv,i)-d(iv,i)*as(iv,i+1))
    307       as(iv,i)=a(iv,i)*x(iv,i)
    308       ds(iv,i)=(e(iv,i)-d(iv,i)*ds(iv,i+1))*x(iv,i)
    309 99998                                              continue
    310 10001 continue
    311                                                    do 99997 iv=1,KDLON
    312       y(iv,1)=ds(iv,1)
    313 99997                                              continue
    314       do 10002 i=2,n
    315                                                    do 99996 iv=1,KDLON
    316       y(iv,i)=ds(iv,i)-as(iv,i)*y(iv,i-1)
    317 99996                                              continue
    318 10002 continue
    319 c.......................................................................
    320 c
    321 c.......................................................................
    322       RETURN
    323       END
     308      do iv=1,KDLON
     309        as(iv,n)=a(iv,n)/b(iv,n)
     310        ds(iv,n)=e(iv,n)/b(iv,n)
     311      enddo
     312      do i=n-1,1,-1
     313        do iv=1,KDLON
     314          x(iv,i)=1.E+0/(b(iv,i)-d(iv,i)*as(iv,i+1))
     315          as(iv,i)=a(iv,i)*x(iv,i)
     316          ds(iv,i)=(e(iv,i)-d(iv,i)*ds(iv,i+1))*x(iv,i)
     317        enddo
     318      enddo
     319      do iv=1,KDLON
     320        y(iv,1)=ds(iv,1)
     321      enddo
     322      do i=2,n
     323        do iv=1,KDLON
     324          y(iv,i)=ds(iv,i)-as(iv,i)*y(iv,i-1)
     325        enddo
     326      enddo
     327c.......................................................................
     328c
     329c.......................................................................
     330
     331      END SUBROUTINE sys3v
     332
     333      END MODULE flusv_mod

  • trunk/LMDZ.MARS/libf/phymars/lwb.F

    r1266 r3757  
     1      module lwb_mod
     2     
     3      implicit none
     4     
     5      contains
     6     
    17      subroutine lwb (kdlon,kflev,tlev,tlay,dt0
    28     .               ,bsurf,btop,blay,blev,dblay,dbsublay)
    … …  
    132138c----------------------------------------------------------------------
    133139
    134       return
    135       end
     140      end subroutine lwb
     141
     142      end module lwb_mod

  • trunk/LMDZ.MARS/libf/phymars/lwdiff.F

    r3726 r3757  
     1      module lwdiff_mod
     2     
     3      implicit none
     4     
     5      contains
     6     
    17      subroutine lwdiff (kdlon,kflev
    28     .         ,pbsur,pbtop,pdbsl
    … …  
    814      use yomlw_h, only: nlaylte
    915      use comcstfi_h, only: pi
     16      use flusv_mod, only: flusv
    1017      IMPLICIT NONE
    1118 
    … …  
    3441C              ---------
    3542C
    36       integer kdlon,kflev
    37       REAL PBSUR(NDLO2,nir), PBTOP(NDLO2,nir)
    38      S  ,  PDBSL(NDLO2,nir,KFLEV*2), PEMIS(NDLO2)
    39 
    40       real PFLUC(NDLO2,2,KFLEV+1)
    41       real tautotal(ndlon,nflev,nir)
    42       real omegtotal(ndlon,nflev,nir), gtotal(ndlon,nflev,nir)
     43      integer,intent(in) :: kdlon,kflev
     44      real,intent(in) :: PBSUR(NDLO2,nir), PBTOP(NDLO2,nir)
     45      real,intent(in) :: PDBSL(NDLO2,nir,KFLEV*2), PEMIS(NDLO2)
     46
     47      real,intent(out) :: PFLUC(NDLO2,2,KFLEV+1)
     48      real,intent(in) :: tautotal(ndlon,nflev,nir)
     49      real,intent(in) :: omegtotal(ndlon,nflev,nir)
     50      real,intent(in) :: gtotal(ndlon,nflev,nir)
    4351
    4452C
    … …  
    6977C                --------------
    7078C
    71  100  CONTINUE
     79
    7280C
    7381C*         1.1     INITIALIZE LAYER CONTRIBUTIONS
    7482C                  ------------------------------
    7583C
    76  110  CONTINUE
     84
    7785C
    7886
    … …  
    8492      enddo
    8593
    86       DO 112 JK = 1 , nlaylte+1
    87         DO 111 JL = 1 , KDLON
     94      DO JK = 1 , nlaylte+1
     95        DO JL = 1 , KDLON
    8896          ZADJD(JL,JK) = 0.
    8997          ZADJU(JL,JK) = 0.
    9098          ZDISD(JL,JK) = 0.
    9199          ZDISU(JL,JK) = 0.
    92  111    CONTINUE
    93  112  CONTINUE
     100        ENDDO
     101      ENDDO
    94102C
    95103C
    … …  
    106114C  ==================================================================
    107115C
    108  200  CONTINUE
     116
    109117C
    110118C     ------------------------------------------------------------------
    … …  
    133141c! boucle sur les  bandes hors CO2
    134142c!
    135       DO 10001 iir=3,nir
     143      DO iir=3,nir
    136144c!
    137145                                            do jl=1,kdlon
    … …  
    202210           enddo
    203211        ENDDO
    204 10001 CONTINUE
     212      ENDDO ! of DO iir=3,nir
    205213
    206214      DO J2=1,nlaylte+1
    … …  
    213221
    214222
    215       END
     223      END SUBROUTINE lwdiff
     224
     225      end module lwdiff_mod

  • trunk/LMDZ.MARS/libf/phymars/lwflux.F

    r3726 r3757  
    1919      use dimradmars_mod, only: ndlo2, nir, ndlon, nuco2, nflev
    2020      use yomlw_h, only: nlaylte, xi, xi_ground, gcp
     21      use lwdiff_mod, only: lwdiff
    2122      implicit none
    2223 

  • trunk/LMDZ.MARS/libf/phymars/lwmain_mod.F

    r3726 r3757  
    2222      use yomlw_h, only: nlaylte, xi
    2323      use lwi_mod, only: lwi
     24      use lwb_mod, only: lwb
     25      use lwu_mod, only: lwu
     26      use lwxd_mod, only: lwxd
     27      use lwxn_mod, only: lwxn
     28      use lwxb_mod, only: lwxb
    2429      use lwflux_mod, only: lwflux
    2530      use callkeys_mod, only: ilwb, ilwd, ilwn

  • trunk/LMDZ.MARS/libf/phymars/lwtt.F

    r1266 r3757  
     1      module lwtt_mod
     2     
     3      implicit none
     4     
     5      contains
     6     
    17      subroutine lwtt (kdlon,u,up,nu,tr)
    28
    … …  
    5965
    6066c----------------------------------------------------------------------
    61       return
    62       end
     67
     68      end subroutine lwtt
     69
     70      end module lwtt_mod

  • trunk/LMDZ.MARS/libf/phymars/lwu.F

    r3726 r3757  
     1      module lwu_mod
     2     
     3      implicit none
     4     
     5      contains
     6     
    17      subroutine lwu (kdlon,kflev
    28     &                ,dp,plev,tlay,aerosol
    … …  
    208214     
    209215c----------------------------------------------------------------------
    210       end
     216      end subroutine lwu
     217
     218      end module lwu_mod

  • trunk/LMDZ.MARS/libf/phymars/lwxb.F

    r1266 r3757  
     1      module lwxb_mod
     2     
     3      implicit none
     4     
     5      contains
     6     
    17      subroutine lwxb (ig0,kdlon,kflev
    28     .                ,emis
    … …  
    3541      use dimradmars_mod, only: ndlo2, nuco2, ndlon, nflev
    3642      use yomlw_h, only: xi, nlaylte
     43      use lwtt_mod, only: lwtt
    3744      implicit none
    3845
    … …  
    212219
    213220c-------------------------------------------------------------------------
    214       return
    215       end
     221      end subroutine lwxb
     222
     223      end module lwxb_mod

  • trunk/LMDZ.MARS/libf/phymars/lwxd.F

    r3726 r3757  
     1      module lwxd_mod
     2     
     3      implicit none
     4     
     5      contains
     6     
    17      subroutine lwxd (ig0,kdlon,kflev,emis
    28     .                ,aer_t,co2_u,co2_up)
    … …  
    3541      use yomlw_h, only: nlaylte, xi, xi_emis
    3642      use callkeys_mod, only: callemis
     43      use lwtt_mod, only: lwtt
    3744      implicit none
    3845
    … …  
    247254
    248255c----------------------------------------------------------------------
    249       end
     256      end subroutine lwxd
     257     
     258      end module lwxd_mod

  • trunk/LMDZ.MARS/libf/phymars/lwxn.F

    r3726 r3757  
     1      module lwxn_mod
     2     
     3      implicit none
     4     
     5      contains
     6     
    17      subroutine lwxn ( ig0,kdlon,kflev
    28     .                , dp
    … …  
    7379      use yomlw_h, only: nlaylte, xi, xi_ground, xi_emis
    7480      use callkeys_mod, only: linear, alphan, ncouche
     81      use lwtt_mod, only: lwtt
    7582      implicit none
    7683
    … …  
    379386      enddo    !  boucle sur jk
    380387c----------------------------------------------------------------------
    381       return
    382       end
    383 
    384 
    385 
     388
     389      end subroutine lwxn
     390
     391      end module lwxn_mod

  • trunk/LMDZ.MARS/libf/phymars/sugwd.F90

    r3754 r3757  
     1MODULE sugwd_mod
     2
     3IMPLICIT NONE
     4
     5CONTAINS
     6
    17      SUBROUTINE SUGWD(nlayer,sigtest)
    28! ==============================================================================
    … …  
    7076      GTSEC=1.E-07    ! Security values for Sub-grid scale anisotropy(OROSETUP,GWSTRESS)
    7177
    72       RETURN
    73       END
     78      END SUBROUTINE SUGWD
     79
     80END MODULE sugwd_mod
Note: See TracChangeset for help on using the changeset viewer.