Ignore:
Timestamp:
Feb 17, 2011, 7:20:44 PM (14 years ago)
Author:
lguez
Message:

Conversion to free source form for "inidissip". "comdissipn.h" is now
valid and equivalent for either free source form or fixed source form.

Added compilation files for Gfortran.

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/inidissip.F90

    r1488 r1489  
    22! $Id$
    33!
    4       SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
    5      *                       tetagdiv,tetagrot,tetatemp             )
    6 c=======================================================================
    7 c   initialisation de la dissipation horizontale
    8 c=======================================================================
    9 c-----------------------------------------------------------------------
    10 c   declarations:
    11 c   -------------
    12 
    13       USE control_mod
    14 
    15       IMPLICIT NONE
    16 #include "dimensions.h"
    17 #include "paramet.h"
    18 #include "comdissipn.h"
    19 #include "comconst.h"
    20 #include "comvert.h"
    21 #include "logic.h"
    22 
    23       LOGICAL lstardis
    24       INTEGER nitergdiv,nitergrot,niterh
    25       REAL    tetagdiv,tetagrot,tetatemp
    26       REAL fact,zvert(llm),zz
    27       REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
    28       REAL ullm,vllm,umin,vmin,zhmin,zhmax
    29       REAL zllm,z1llm
    30 
    31       INTEGER l,ij,idum,ii
    32       REAL tetamin
    33       REAL pseudoz
    34 
    35       REAL ran1
    36 
    37 
    38 c-----------------------------------------------------------------------
    39 c
    40 c   calcul des valeurs propres des operateurs par methode iterrative:
    41 c   -----------------------------------------------------------------
    42 
    43       crot     = -1.
    44       cdivu    = -1.
    45       cdivh    = -1.
    46 
    47 c   calcul de la valeur propre de divgrad:
    48 c   --------------------------------------
    49       idum = 0
    50       DO l = 1, llm
    51        DO ij = 1, ip1jmp1
     4SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
     5     tetagdiv,tetagrot,tetatemp             )
     6  !=======================================================================
     7  !   initialisation de la dissipation horizontale
     8  !=======================================================================
     9  !-----------------------------------------------------------------------
     10  !   declarations:
     11  !   -------------
     12
     13  USE control_mod
     14
     15  IMPLICIT NONE
     16  include "dimensions.h"
     17  include "paramet.h"
     18  include "comdissipn.h"
     19  include "comconst.h"
     20  include "comvert.h"
     21  include "logic.h"
     22
     23  LOGICAL lstardis
     24  INTEGER nitergdiv,nitergrot,niterh
     25  REAL    tetagdiv,tetagrot,tetatemp
     26  REAL fact,zvert(llm),zz
     27  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
     28  REAL ullm,vllm,umin,vmin,zhmin,zhmax
     29  REAL zllm,z1llm
     30
     31  INTEGER l,ij,idum,ii
     32  REAL tetamin
     33  REAL pseudoz
     34
     35  REAL ran1
     36
     37
     38  !-----------------------------------------------------------------------
     39  !
     40  !   calcul des valeurs propres des operateurs par methode iterrative:
     41  !   -----------------------------------------------------------------
     42
     43  crot     = -1.
     44  cdivu    = -1.
     45  cdivh    = -1.
     46
     47  !   calcul de la valeur propre de divgrad:
     48  !   --------------------------------------
     49  idum = 0
     50  DO l = 1, llm
     51     DO ij = 1, ip1jmp1
    5252        deltap(ij,l) = 1.
    53        ENDDO
    54       ENDDO
    55 
    56       idum  = -1
    57       zh(1) = RAN1(idum)-.5
    58       idum  = 0
    59       DO ij = 2, ip1jmp1
    60         zh(ij) = RAN1(idum) -.5
    61       ENDDO
    62 
    63       CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
    64 
    65       CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    66 
    67       IF ( zhmin .GE. zhmax  )     THEN
    68          PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
    69          STOP'probleme generateur alleatoire dans inidissip'
    70       ENDIF
    71 
    72       zllm = ABS( zhmax )
    73       DO l = 1,50
    74          IF(lstardis) THEN
    75             CALL divgrad2(1,zh,deltap,niterh,zh)
    76          ELSE
    77             CALL divgrad (1,zh,niterh,zh)
    78          ENDIF
    79 
    80         CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    81 
    82          zllm  = ABS( zhmax )
    83          z1llm = 1./zllm
    84          DO ij = 1,ip1jmp1
    85             zh(ij) = zh(ij)* z1llm
    86          ENDDO
    87       ENDDO
    88 
    89       IF(lstardis) THEN
    90          cdivh = 1./ zllm
    91       ELSE
    92          cdivh = zllm ** ( -1./niterh )
    93       ENDIF
    94 
    95 c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
    96 c   -----------------------------------------------------------------
    97       print*,'calcul des valeurs propres'
    98 
    99       DO  20  ii = 1, 2
    100 c
    101          DO ij = 1, ip1jmp1
    102            zu(ij)  = RAN1(idum) -.5
    103          ENDDO
    104          CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
    105          DO ij = 1, ip1jm
    106             zv(ij) = RAN1(idum) -.5
    107          ENDDO
    108          CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
    109 
    110          CALL minmax(iip1*jjp1,zu,umin,ullm )
    111          CALL minmax(iip1*jjm, zv,vmin,vllm )
    112 
    113          ullm = ABS ( ullm )
    114          vllm = ABS ( vllm )
    115 
    116          DO  5  l = 1, 50
    117             IF(ii.EQ.1) THEN
    118 ccccc             CALL covcont( 1,zu,zv,zu,zv )
    119                IF(lstardis) THEN
    120                   CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
    121                ELSE
    122                   CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
    123                ENDIF
    124             ELSE
    125                IF(lstardis) THEN
    126                   CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
    127                ELSE
    128                   CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
    129                ENDIF
    130             ENDIF
    131 
    132             CALL minmax(iip1*jjp1,zu,umin,ullm )
    133             CALL minmax(iip1*jjm, zv,vmin,vllm )
    134 
    135             ullm = ABS  ( ullm )
    136             vllm = ABS  ( vllm )
    137 
    138             zllm  = MAX( ullm,vllm )
    139             z1llm = 1./ zllm
    140             DO ij = 1, ip1jmp1
    141               zu(ij) = zu(ij)* z1llm
    142             ENDDO
    143             DO ij = 1, ip1jm
    144                zv(ij) = zv(ij)* z1llm
    145             ENDDO
    146  5       CONTINUE
    147 
    148          IF ( ii.EQ.1 ) THEN
    149             IF(lstardis) THEN
    150                cdivu  = 1./zllm
    151             ELSE
    152                cdivu  = zllm **( -1./nitergdiv )
    153             ENDIF
    154          ELSE
    155             IF(lstardis) THEN
    156                crot   = 1./ zllm
    157             ELSE
    158                crot   = zllm **( -1./nitergrot )
    159             ENDIF
    160          ENDIF
    161 
    162  20   CONTINUE
    163 
    164 c   petit test pour les operateurs non star:
    165 c   ----------------------------------------
    166 
    167 c     IF(.NOT.lstardis) THEN
    168          fact    = rad*24./REAL(jjm)
    169          fact    = fact*fact
    170          PRINT*,'coef u ', fact/cdivu, 1./cdivu
    171          PRINT*,'coef r ', fact/crot , 1./crot
    172          PRINT*,'coef h ', fact/cdivh, 1./cdivh
    173 c     ENDIF
    174 
    175 c-----------------------------------------------------------------------
    176 c   variation verticale du coefficient de dissipation:
    177 c   --------------------------------------------------
    178 
    179       if (ok_strato .and. llm==39) then
    180          do l=1,llm
    181             pseudoz=8.*log(preff/presnivs(l))
    182             zvert(l)=1+
    183      s      (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2.
    184      s      *(dissip_factz-1.)
    185          enddo
    186       else
    187          DO l=1,llm
    188             zvert(l)=1.
    189          ENDDO
    190          fact=2.
    191          DO l = 1, llm
    192             zz      = 1. - preff/presnivs(l)
    193             zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
    194          ENDDO
    195       endif
    196 
    197 
    198       PRINT*,'Constantes de temps de la diffusion horizontale'
    199 
    200       tetamin =  1.e+6
    201 
    202       DO l=1,llm
    203         tetaudiv(l)   = zvert(l)/tetagdiv
    204         tetaurot(l)   = zvert(l)/tetagrot
    205         tetah(l)      = zvert(l)/tetatemp
    206 
    207         IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
    208         IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
    209         IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
    210       ENDDO
    211 
    212       PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
    213       idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
    214       PRINT *,' INIDI tetamin idissip ',tetamin,idissip
    215       idissip = MAX(iperiod,idissip)
    216       dtdiss  = idissip * dtvr
    217       PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
    218 
    219       DO l = 1,llm
    220          PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
    221      *                   dtdiss*tetah(l)
    222       ENDDO
    223 
    224 c
    225       RETURN
    226       END
     53     ENDDO
     54  ENDDO
     55
     56  idum  = -1
     57  zh(1) = RAN1(idum)-.5
     58  idum  = 0
     59  DO ij = 2, ip1jmp1
     60     zh(ij) = RAN1(idum) -.5
     61  ENDDO
     62
     63  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
     64
     65  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
     66
     67  IF ( zhmin .GE. zhmax  )     THEN
     68     PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
     69     STOP'probleme generateur alleatoire dans inidissip'
     70  ENDIF
     71
     72  zllm = ABS( zhmax )
     73  DO l = 1,50
     74     IF(lstardis) THEN
     75        CALL divgrad2(1,zh,deltap,niterh,zh)
     76     ELSE
     77        CALL divgrad (1,zh,niterh,zh)
     78     ENDIF
     79
     80     CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
     81
     82     zllm  = ABS( zhmax )
     83     z1llm = 1./zllm
     84     DO ij = 1,ip1jmp1
     85        zh(ij) = zh(ij)* z1llm
     86     ENDDO
     87  ENDDO
     88
     89  IF(lstardis) THEN
     90     cdivh = 1./ zllm
     91  ELSE
     92     cdivh = zllm ** ( -1./niterh )
     93  ENDIF
     94
     95  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
     96  !   -----------------------------------------------------------------
     97  print*,'calcul des valeurs propres'
     98
     99  DO    ii = 1, 2
     100     !
     101     DO ij = 1, ip1jmp1
     102        zu(ij)  = RAN1(idum) -.5
     103     ENDDO
     104     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
     105     DO ij = 1, ip1jm
     106        zv(ij) = RAN1(idum) -.5
     107     ENDDO
     108     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
     109
     110     CALL minmax(iip1*jjp1,zu,umin,ullm )
     111     CALL minmax(iip1*jjm, zv,vmin,vllm )
     112
     113     ullm = ABS ( ullm )
     114     vllm = ABS ( vllm )
     115
     116     DO    l = 1, 50
     117        IF(ii.EQ.1) THEN
     118           !cccc             CALL covcont( 1,zu,zv,zu,zv )
     119           IF(lstardis) THEN
     120              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
     121           ELSE
     122              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
     123           ENDIF
     124        ELSE
     125           IF(lstardis) THEN
     126              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
     127           ELSE
     128              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
     129           ENDIF
     130        ENDIF
     131
     132        CALL minmax(iip1*jjp1,zu,umin,ullm )
     133        CALL minmax(iip1*jjm, zv,vmin,vllm )
     134
     135        ullm = ABS  ( ullm )
     136        vllm = ABS  ( vllm )
     137
     138        zllm  = MAX( ullm,vllm )
     139        z1llm = 1./ zllm
     140        DO ij = 1, ip1jmp1
     141           zu(ij) = zu(ij)* z1llm
     142        ENDDO
     143        DO ij = 1, ip1jm
     144           zv(ij) = zv(ij)* z1llm
     145        ENDDO
     146     end DO
     147
     148     IF ( ii.EQ.1 ) THEN
     149        IF(lstardis) THEN
     150           cdivu  = 1./zllm
     151        ELSE
     152           cdivu  = zllm **( -1./nitergdiv )
     153        ENDIF
     154     ELSE
     155        IF(lstardis) THEN
     156           crot   = 1./ zllm
     157        ELSE
     158           crot   = zllm **( -1./nitergrot )
     159        ENDIF
     160     ENDIF
     161
     162  end DO
     163
     164  !   petit test pour les operateurs non star:
     165  !   ----------------------------------------
     166
     167  !     IF(.NOT.lstardis) THEN
     168  fact    = rad*24./REAL(jjm)
     169  fact    = fact*fact
     170  PRINT*,'coef u ', fact/cdivu, 1./cdivu
     171  PRINT*,'coef r ', fact/crot , 1./crot
     172  PRINT*,'coef h ', fact/cdivh, 1./cdivh
     173  !     ENDIF
     174
     175  !-----------------------------------------------------------------------
     176  !   variation verticale du coefficient de dissipation:
     177  !   --------------------------------------------------
     178
     179  if (ok_strato .and. llm==39) then
     180     do l=1,llm
     181        pseudoz=8.*log(preff/presnivs(l))
     182        zvert(l)=1+ &
     183             (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
     184             *(dissip_factz-1.)
     185     enddo
     186  else
     187     DO l=1,llm
     188        zvert(l)=1.
     189     ENDDO
     190     fact=2.
     191     DO l = 1, llm
     192        zz      = 1. - preff/presnivs(l)
     193        zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
     194     ENDDO
     195  endif
     196
     197
     198  PRINT*,'Constantes de temps de la diffusion horizontale'
     199
     200  tetamin =  1.e+6
     201
     202  DO l=1,llm
     203     tetaudiv(l)   = zvert(l)/tetagdiv
     204     tetaurot(l)   = zvert(l)/tetagrot
     205     tetah(l)      = zvert(l)/tetatemp
     206
     207     IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
     208     IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
     209     IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
     210  ENDDO
     211
     212  PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
     213  idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
     214  PRINT *,' INIDI tetamin idissip ',tetamin,idissip
     215  idissip = MAX(iperiod,idissip)
     216  dtdiss  = idissip * dtvr
     217  PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
     218
     219  DO l = 1,llm
     220     PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
     221          dtdiss*tetah(l)
     222  ENDDO
     223
     224END SUBROUTINE inidissip
Note: See TracChangeset for help on using the changeset viewer.