Ignore:
Timestamp:
Sep 19, 2025, 2:53:16 PM (4 months ago)
Author:
jbclement
Message:

Mars PCM 1D:

  • "profile.F" is transformed into a module "profile_temp_mod.F90" with F90 format. The file "profile" where the temperature is read is renamed into "profile_temp" to avoid any confusion.
  • Small corrections and improvements throughout the code and scripts in the deftank.

JBC

File:
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/profile_temp_mod.F90

    r3917 r3918  
    1 !      SUBROUTINE profile(unit,nlev,zkm,temp)
    2       SUBROUTINE profile(nlev,zkm,temp)
    3 ! to use  'getin'
    4       USE ioipsl_getincom
    5       IMPLICIT NONE
    6 c=======================================================================
    7 c     Subroutine utilisee dans le modele 1-D  "testphys1d"
    8 c     pour l'initialisation du profil atmospherique
    9 c=======================================================================
    10 c
    11 c   differents profils d'atmospheres. T=f(z)
    12 c   entree:
    13 c     nlev    nombre de niveaux
    14 c     zkm     alititudes en km
    15 c     ichoice choix de l'atmosphere:
    16 c             1 Temperature constante
    17 c             2 profil Savidjari
    18 c             3 Lindner (profil polaire)
    19 c             4 Inversion pour Francois
    20 c             5 Seiff (moyen)
    21 c             6 T constante + perturbation gauss (level) (christophe 10/98)
    22 c             7 T constante + perturbation gauss   (km)  (christophe 10/98)
    23 c             8 Lecture du profile dans un fichier ASCII (profile)
    24 c     tref    temperature de reference
    25 c     isin    ajout d'une perturbation (isin=1)
    26 c     pic     pic perturbation gauss pour ichoice = 6 ou 7
    27 c     largeur largeur de la perturbation gauss pour ichoice = 6 ou 7
    28 c     hauteur hauteur de la perturbation gauss pour ichoice = 6 ou 7
    29 c
    30 c   sortie:
    31 c     temp    temperatures en K
    32 c     
    33 c=======================================================================
    34 c-----------------------------------------------------------------------
    35 c   declarations:
    36 c   -------------
    37 
    38 c   arguments:
    39 c   ----------
    40 
    41        INTEGER nlev
    42        REAL zkm(nlev),temp(nlev)
    43 
    44 c   local:
    45 c   ------
    46 
    47       INTEGER il,ichoice,nseiff,iseiff,isin,iter
    48       REAL pi
    49       PARAMETER(nseiff=37)
    50       REAL tref,t1,t2,t3,ww
    51       REAL tseiff(nseiff)
    52       DATA tseiff/214.,213.8,213.4,212.4,209.3,205.,201.4,197.8,
    53      $           194.6,191.4,188.2,185.2,182.5,180.,177.5,175,
    54      $           172.5,170.,167.5,164.8,162.4,160.,158.,156.,
    55      $           154.1,152.2,150.3,148.7,147.2,145.7,144.2,143.,
    56      $           142.,141.,140,139.5,139./
    57       REAL pic,largeur
    58       REAL hauteur,tmp
    59 
    60 c-----------------------------------------------------------------------
    61 c   read input profile type:
    62 c--------------------------
    63 
    64       ichoice=1 ! default value for ichoice
    65       call getin("ichoice",ichoice)
    66       tref=200 ! default value for tref
    67       call getin("tref",tref)
    68       isin=0 ! default value for isin (=0 means no perturbation)
    69       call getin("isin",isin)
    70       pic=26.522 ! default value for pic
    71       call getin("pic",pic)
    72       largeur=10 ! default value for largeur
    73       call getin("largeur",largeur)
    74       hauteur=30 ! default value for hauteur
    75       call getin("hauteur",hauteur)
    76 
    77 c-----------------------------------------------------------------------
    78 c   ichoice=1 temperature constante:
    79 c   --------------------------------
    80 
    81       IF(ichoice.EQ.1) THEN
    82          DO il=1,nlev
    83             temp(il)=tref
    84          ENDDO
    85 
    86 c-----------------------------------------------------------------------
    87 c   ichoice=2 savidjari:
    88 c   --------------------
    89 
    90       ELSE IF(ichoice.EQ.2) THEN
    91          DO il=1,nlev
    92             temp(il)=AMAX1(219.-2.5*zkm(il),140.)
    93          ENDDO
    94 
    95 c-----------------------------------------------------------------------
    96 c   ichoice=3 Lindner:
    97 c   ------------------
    98 
    99       ELSE IF(ichoice.EQ.3) THEN
    100          DO il=1,nlev
    101             IF(zkm(il).LT.2.5) THEN
    102                temp(il)=150.+30.*zkm(il)/2.5
    103             ELSE IF(zkm(il).LT.5.) THEN
    104                temp(il)=180.
    105             ELSE
    106                temp(il)=AMAX1(180.-2.*(zkm(il)-5.),130.)
    107             ENDIF
    108          ENDDO
    109 
    110 c-----------------------------------------------------------------------
    111 c   ichoice=4 Inversion pour Francois:
    112 c   ------------------
    113 
    114       ELSE IF(ichoice.EQ.4) THEN
    115          DO il=1,nlev
    116             IF(zkm(il).LT.20.) THEN
    117                temp(il)=135.
    118             ELSE
    119                temp(il)=AMIN1(135.+5.*(zkm(il)-20.),200.)
    120             ENDIF
    121          ENDDO
    122 
    123 
    124 c-----------------------------------------------------------------------
    125 c   ichoice=5 Seiff:
    126 c   ----------------
    127 
    128       ELSE IF(ichoice.EQ.5) THEN
    129          DO il=1,nlev
    130             iseiff=INT(zkm(il)/2.)+1
    131             IF(iseiff.LT.nseiff-1) THEN
    132                temp(il)=tseiff(iseiff)+(zkm(il)-2.*(iseiff-1))*
    133      $         (tseiff(iseiff+1)-tseiff(iseiff))/2.
    134             ELSE
    135                temp(il)=tseiff(nseiff)
    136             ENDIF
    137          ENDDO
    138 c IF(ichoice.EQ.6) THEN
    139 c           DO iter=1,6
    140 c           t2=temp(1)
    141 c           t3=temp(2)
    142 c           DO il=2,nlev-1
    143 c              t1=t2
    144 c              t2=t3
    145 c              t3=temp(il+1)
    146 c              ww=tanh(zkm(il)/20.)
    147 c              ww=ww*ww*ww
    148 c              temp(il)=t2+ww*.5*(t1-2.*t2+t3)
    149 c           ENDDO
    150 c           ENDDO
    151 c        ENDIF
    152 
    153 c-----------------------------------------------------------------------
    154 c   ichoice=6
    155 c   ---------
    156 
    157       ELSE IF(ichoice.EQ.6) THEN
    158       DO il=1,nlev
     1MODULE profile_temp_mod
     2
     3implicit none
     4
     5CONTAINS
     6
     7SUBROUTINE profile_temp(nlev,zkm,temp)
     8
     9use ioipsl_getincom, only: getin ! To use 'getin'
     10use comcstfi_h,      only: pi
     11     
     12IMPLICIT NONE
     13
     14!=======================================================================
     15!     Subroutine utilisee dans le modele 1-D "testphys1d"
     16!     pour l'initialisation du profil atmospherique
     17!=======================================================================
     18!
     19!   differents profils d'atmospheres. T=f(z)
     20!   entree:
     21!     nlev    nombre de niveaux
     22!     zkm     alititudes en km
     23!     ichoice choix de l'atmosphere:
     24!             1 Temperature constante
     25!             2 profil Savidjari
     26!             3 Lindner (profil polaire)
     27!             4 Inversion pour Francois
     28!             5 Seiff (moyen)
     29!             6 T constante + perturbation gauss (level) (christophe 10/98)
     30!             7 T constante + perturbation gauss    (km) (christophe 10/98)
     31!             8 Lecture du profile dans un fichier ASCII (profile_temp)
     32!     tref    temperature de reference
     33!     isin    ajout d'une perturbation (isin=1)
     34!     pic     pic perturbation gauss pour ichoice = 6 ou 7
     35!     largeur largeur de la perturbation gauss pour ichoice = 6 ou 7
     36!     hauteur hauteur de la perturbation gauss pour ichoice = 6 ou 7
     37!
     38!   sortie:
     39!     temp    temperatures en K
     40!     
     41!=======================================================================
     42!-----------------------------------------------------------------------
     43!   declarations:
     44!   -------------
     45
     46!   arguments:
     47!   ----------
     48
     49INTEGER,               intent(in) :: nlev
     50REAL, dimension(nlev), intent(in) :: zkm
     51REAL, dimension(nlev), intent(out) :: temp
     52
     53!   local:
     54!   ------
     55
     56INTEGER :: il,ichoice,iseiff,isin,iter,ierr
     57INTEGER, PARAMETER :: nseiff = 37
     58REAL :: tref,t1,t2,t3,ww
     59REAL, dimension(nseiff) :: tseiff
     60DATA tseiff/214. ,213.8,213.4,212.4,209.3,205. ,201.4,197.8, &
     61            194.6,191.4,188.2,185.2,182.5,180. ,177.5,175  , &
     62            172.5,170. ,167.5,164.8,162.4,160. ,158. ,156. , &
     63            154.1,152.2,150.3,148.7,147.2,145.7,144.2,143. , &
     64            142. ,141. ,140  ,139.5,139./
     65REAL :: pic,largeur
     66REAL :: hauteur,tmp
     67
     68!-----------------------------------------------------------------------
     69!   read input profile type:
     70!--------------------------
     71
     72ichoice=1 ! default value for ichoice
     73call getin("ichoice",ichoice)
     74tref=200 ! default value for tref
     75call getin("tref",tref)
     76isin=0 ! default value for isin (=0 means no perturbation)
     77call getin("isin",isin)
     78pic=26.522 ! default value for pic
     79call getin("pic",pic)
     80largeur=10 ! default value for largeur
     81call getin("largeur",largeur)
     82hauteur=30 ! default value for hauteur
     83call getin("hauteur",hauteur)
     84
     85!-----------------------------------------------------------------------
     86!   ichoice=1 temperature constante:
     87!   --------------------------------
     88
     89IF(ichoice.EQ.1) THEN
     90    DO il=1,nlev
     91        temp(il)=tref
     92    ENDDO
     93
     94!-----------------------------------------------------------------------
     95!   ichoice=2 savidjari:
     96!   --------------------
     97
     98ELSE IF(ichoice.EQ.2) THEN
     99    DO il=1,nlev
     100        temp(il)=AMAX1(219.-2.5*zkm(il),140.)
     101    ENDDO
     102
     103!-----------------------------------------------------------------------
     104!   ichoice=3 Lindner:
     105!   ------------------
     106
     107ELSE IF(ichoice.EQ.3) THEN
     108    DO il=1,nlev
     109        IF(zkm(il).LT.2.5) THEN
     110            temp(il)=150.+30.*zkm(il)/2.5
     111        ELSE IF(zkm(il).LT.5.) THEN
     112            temp(il)=180.
     113        ELSE
     114            temp(il)=AMAX1(180.-2.*(zkm(il)-5.),130.)
     115        ENDIF
     116    ENDDO
     117
     118!-----------------------------------------------------------------------
     119!   ichoice=4 Inversion pour Francois:
     120!   ------------------
     121
     122ELSE IF(ichoice.EQ.4) THEN
     123    DO il=1,nlev
     124        IF(zkm(il).LT.20.) THEN
     125            temp(il)=135.
     126        ELSE
     127            temp(il)=AMIN1(135.+5.*(zkm(il)-20.),200.)
     128        ENDIF
     129    ENDDO
     130
     131!-----------------------------------------------------------------------
     132!   ichoice=5 Seiff:
     133!   ----------------
     134
     135ELSE IF(ichoice.EQ.5) THEN
     136    DO il=1,nlev
     137        iseiff=INT(zkm(il)/2.)+1
     138        IF(iseiff.LT.nseiff-1) THEN
     139            temp(il)=tseiff(iseiff)+(zkm(il)-2.*(iseiff-1))*(tseiff(iseiff+1)-tseiff(iseiff))/2.
     140        ELSE
     141            temp(il)=tseiff(nseiff)
     142        ENDIF
     143    ENDDO
     144! IF(ichoice.EQ.6) THEN
     145!    DO iter=1,6
     146!        t2=temp(1)
     147!        t3=temp(2)
     148!        DO il=2,nlev-1
     149!            t1=t2
     150!            t2=t3
     151!            t3=temp(il+1)
     152!            ww=tanh(zkm(il)/20.)
     153!            ww=ww*ww*ww
     154!            temp(il)=t2+ww*.5*(t1-2.*t2+t3)
     155!        ENDDO
     156!    ENDDO
     157!ENDIF
     158
     159!-----------------------------------------------------------------------
     160!   ichoice=6
     161!   ---------
     162
     163ELSE IF(ichoice.EQ.6) THEN
     164    DO il=1,nlev
    159165        tmp=il-pic
    160166        temp(il)=tref + hauteur*exp(-tmp*tmp/largeur/largeur)
    161       ENDDO
    162 
    163 
    164 c-----------------------------------------------------------------------
    165 c   ichoice=7
    166 c   ---------
    167 
    168       ELSE IF(ichoice.EQ.7) THEN
    169       DO il=1,nlev
     167    ENDDO
     168
     169!-----------------------------------------------------------------------
     170!   ichoice=7
     171!   ---------
     172
     173ELSE IF(ichoice.EQ.7) THEN
     174    DO il=1,nlev
    170175        tmp=zkm(il)-pic
    171176        temp(il)=tref + hauteur*exp(-tmp*tmp*4/largeur/largeur)
    172       ENDDO
    173 
    174 c-----------------------------------------------------------------------
    175 c   ichoice=8
    176 c   ---------
    177 
    178       ELSE IF(ichoice.GE.8) THEN
    179       OPEN(11,file='profile',status='old',form='formatted',err=101)
    180       DO il=1,nlev
    181         READ (11,*) temp(il)
    182       ENDDO
    183 
    184       GOTO 201
    185 101   STOP'fichier profile inexistant'
    186 201   CONTINUE
    187       CLOSE(10)
    188 
    189 c-----------------------------------------------------------------------
    190 
    191       ENDIF
    192 
    193 c-----------------------------------------------------------------------
    194 c   rajout eventuel d'une perturbation:
    195 c   -----------------------------------
    196 
    197       IF(isin.EQ.1) THEN
    198          pi=2.*ASIN(1.)
    199          DO il=1,nlev
    200 c       if (nlev.EQ.501) then
    201 c         if (zkm(il).LE.70.5) then
    202 c       temp(il)=temp(il)+(1.-1000./(1000+zkm(il)*zkm(il)))*(
    203 c    s      6.*SIN(zkm(il)*pi/6.)+9.*SIN(zkm(il)*pi/10.3) )
    204 c         endif
    205 c       else
    206         temp(il)=temp(il)+(1.-1000./(1000+zkm(il)*zkm(il)))*(
    207      s      6.*SIN(zkm(il)*pi/6.)+9.*SIN(zkm(il)*pi/10.3) )
    208 c       endif
    209          ENDDO
    210       ENDIF
    211 
    212 
    213 c-----------------------------------------------------------------------
    214 c   Ecriture du profil de temperature dans un fichier profile.out
    215 c   -------------------------------------------------------------
    216 
    217 
    218       OPEN(12,file='profile.out',form='formatted')
    219           DO il=1,nlev
    220             write(12,*) temp(il)
    221           ENDDO
    222       CLOSE(12)
    223 
    224       RETURN
    225       END
     177    ENDDO
     178
     179!-----------------------------------------------------------------------
     180!   ichoice=8
     181!   ---------
     182
     183ELSE IF(ichoice.GE.8) THEN
     184    OPEN(11,file='profile_temp',status='old',form='formatted',iostat=ierr)
     185    if (ierr == 0) then
     186        DO il=1,nlev
     187            READ (11,*) temp(il)
     188        ENDDO
     189    else
     190        error stop 'File "profile_temp" not found!'
     191    endif
     192    CLOSE(11)
     193
     194!-----------------------------------------------------------------------
     195ENDIF
     196
     197!-----------------------------------------------------------------------
     198!   rajout eventuel d'une perturbation:
     199!   -----------------------------------
     200
     201IF (isin.EQ.1) THEN
     202    DO il=1,nlev
     203!        if (nlev.EQ.501) then
     204!            if (zkm(il).LE.70.5) then
     205!                temp(il)=temp(il)+(1.-1000./(1000+zkm(il)*zkm(il)))*( 6.*SIN(zkm(il)*pi/6.)+9.*SIN(zkm(il)*pi/10.3) )
     206!            endif
     207!        else
     208        temp(il)=temp(il)+(1.-1000./(1000+zkm(il)*zkm(il)))*( 6.*SIN(zkm(il)*pi/6.)+9.*SIN(zkm(il)*pi/10.3) )
     209!        endif
     210    ENDDO
     211ENDIF
     212
     213!-----------------------------------------------------------------------
     214!   Ecriture du profil de temperature dans un fichier profile_temp.out
     215!   ------------------------------------------------------------------
     216
     217OPEN(12,file='profile_temp.out',status='unknown',form='formatted')
     218    DO il=1,nlev
     219        write(12,*) temp(il)
     220    ENDDO
     221CLOSE(12)
     222
     223END SUBROUTINE profile_temp
     224
     225END MODULE profile_temp_mod
Note: See TracChangeset for help on using the changeset viewer.