Changeset 5272 for LMDZ6/trunk/libf/phylmd/dyn1d

Timestamp:

Oct 24, 2024, 5:53:15 PM (4 months ago)

Author:

abarral

Message:

Turn paramet.h into a module

Location:

LMDZ6/trunk/libf/phylmd/dyn1d

Files:

• 1DUTILS.h (modified) (47 diffs)
• lmdz1d.F90 (modified) (1 diff)
• paramet_mod_h.f90 (moved) (moved from LMDZ6/trunk/libf/phylmd/dyn1d/paramet.h) (1 diff)

LMDZ6/trunk/libf/phylmd/dyn1d/1DUTILS.h

@@ -1 @@
-INCLUDE "conf_gcm.f90"
+#include "conf_gcm.f90"
 
 !
@@ -18 +18 @@
 !   --------------
 
-INCLUDE "compar1d.h"
-INCLUDE "flux_arp.h"
-INCLUDE "tsoilnudge.h"
-INCLUDE "fcg_gcssold.h"
-INCLUDE "fcg_racmo.h"
+#include "compar1d.h"
+#include "flux_arp.h"
+#include "tsoilnudge.h"
+#include "fcg_gcssold.h"
+#include "fcg_racmo.h"
 !
 !
@@ -29 +29 @@
 
 !      CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
-      
+
 !
 !  -------------------------------------------------------------------
@@ -42 +42 @@
 !Config  Desc = unite de fichier pour les impressions
 !Config  Def  = 6
-!Config  Help = unite de fichier pour les impressions 
+!Config  Help = unite de fichier pour les impressions
 !Config         (defaut sortie standard = 6)
       lunout=6
@@ -74 +74 @@
 !!Config  Help = 0 ==> forcing_les = .true.
 !             initial profiles from file prof.inp.001
-!             no forcing by LS convergence ; 
+!             no forcing by LS convergence ;
 !             surface temperature imposed ;
 !             radiative cooling may be imposed (iflag_radia=0 in physiq.def)
 !         = 1 ==> forcing_radconv = .true.
-!             idem forcing_type = 0, but the imposed radiative cooling 
-!             is set to 0 (hence, if iflag_radia=0 in physiq.def, 
+!             idem forcing_type = 0, but the imposed radiative cooling
+!             is set to 0 (hence, if iflag_radia=0 in physiq.def,
 !             then there is no radiative cooling at all)
 !         = 2 ==> forcing_toga = .true.
-!             initial profiles from TOGA-COARE IFA files 
-!             LS convergence and SST imposed from TOGA-COARE IFA files 
+!             initial profiles from TOGA-COARE IFA files
+!             LS convergence and SST imposed from TOGA-COARE IFA files
 !         = 3 ==> forcing_GCM2SCM = .true.
 !             initial profiles from the GCM output
 !             LS convergence imposed from the GCM output
 !         = 4 ==> forcing_twpi = .true.
-!             initial profiles from TWPICE nc files 
-!             LS convergence and SST imposed from TWPICE nc files 
+!             initial profiles from TWPICE nc files
+!             LS convergence and SST imposed from TWPICE nc files
 !         = 5 ==> forcing_rico = .true.
 !             initial profiles from RICO idealized
-!             LS convergence imposed from  RICO (cst) 
+!             LS convergence imposed from  RICO (cst)
 !         = 6 ==> forcing_amma = .true.
 !         = 10 ==> forcing_case = .true.
-!             initial profiles from case.nc file 
+!             initial profiles from case.nc file
 !         = 40 ==> forcing_GCSSold = .true.
 !             initial profile from GCSS file
@@ -105 +105 @@
 !             Radiation has to be computed interactively
 !         = 60 ==> forcing_astex = .true.
-!             initial profiles from file: see prof.inp.001 
+!             initial profiles from file: see prof.inp.001
 !             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
 !             Radiation has to be computed interactively
 !         = 61 ==> forcing_armcu = .true.
-!             initial profiles from file: see prof.inp.001 
+!             initial profiles from file: see prof.inp.001
 !             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
 !             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
-!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 
+!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s
 !             Radiation to be switched off
 !         > 100 ==> forcing_case = .true. or forcing_case2 = .true.
-!             initial profiles from case.nc file 
+!             initial profiles from case.nc file
 !
        forcing_type = 0
        CALL getin('forcing_type',forcing_type)
          imp_fcg_gcssold   = .false.
-         ts_fcg_gcssold    = .false.  
-         Tp_fcg_gcssold    = .false.  
-         Tp_ini_gcssold    = .false.  
-         xTurb_fcg_gcssold = .false. 
+         ts_fcg_gcssold    = .false.
+         Tp_fcg_gcssold    = .false.
+         Tp_ini_gcssold    = .false.
+         xTurb_fcg_gcssold = .false.
         IF (forcing_type .eq.40) THEN
           CALL getin('imp_fcg',imp_fcg_gcssold)
@@ -261 +261 @@
 !Config  Desc = meaningless in this  case
 !Config  Def  = 0.
-!Config  Help = 
+!Config  Help =
        time_ini = 0.
        CALL getin('time_ini',time_ini)
@@ -277 +277 @@
 !Config  Desc = Grid cell area
 !Config  Def  = 1.e11
-!Config  Help = 
+!Config  Help =
        airefi = 1.e11
        CALL getin('airephy',airefi)
@@ -298 +298 @@
 !Config  Desc = surface pressure
 !Config  Def  = 102400.
-!Config  Help = 
+!Config  Help =
        psurf = 102400.
        CALL getin('psurf',psurf)
@@ -305 +305 @@
 !Config  Desc = surface altitude
 !Config  Def  = 0.
-!Config  Help = 
+!Config  Help =
        zsurf = 0.
        CALL getin('zsurf',zsurf)
-! EV pour accord avec format standard       
+! EV pour accord avec format standard
        CALL getin('zorog',zsurf)
 
@@ -340 +340 @@
 !Config  Desc = ???
 !Config  Def  = 0.0 0.0
-!Config  Help = 
+!Config  Help =
        wtsurf = 0.0
        wqsurf = 0.0
@@ -349 +349 @@
 !Config  Desc = albedo
 !Config  Def  = 0.09
-!Config  Help = 
+!Config  Help =
        albedo = 0.09
        CALL getin('albedo',albedo)
@@ -356 +356 @@
 !Config  Desc = age de la neige
 !Config  Def  = 30.0
-!Config  Help = 
+!Config  Help =
        xagesno = 30.0
        CALL getin('agesno',xagesno)
@@ -363 +363 @@
 !Config  Desc = age de la neige
 !Config  Def  = 30.0
-!Config  Help = 
+!Config  Help =
        restart_runoff = 0.0
        CALL getin('restart_runoff',restart_runoff)
@@ -370 +370 @@
 !Config  Desc = initial bucket water content (kg/m2) when land (5std)
 !Config  Def  = 30.0
-!Config  Help = 
+!Config  Help =
        qsolinp = 1.
        CALL getin('qsolinp',qsolinp)
@@ -379 +379 @@
 !Config  Desc = beta for actual evaporation when prescribed
 !Config  Def  = 1.0
-!Config  Help = 
+!Config  Help =
        betaevap = 1.
-       CALL getin('betaevap',betaevap)      
+       CALL getin('betaevap',betaevap)
 
 !Config  Key  = zpicinp
@@ -689 +689 @@
       real :: q(klon,klev,nqtot),omega2(klon,klev)
 !      real :: ug(klev),vg(klev),fcoriolis
-      real :: phis(klon) 
+      real :: phis(klon)
 
 !   Variables locales pour NetCDF:
@@ -719 +719 @@
 !
       CALL get_var("controle",tab_cntrl)
-       
+
 
       im         = tab_cntrl(1)
@@ -755 +755 @@
         fxyhypb = .false.
         ysinus  = .false.
-        IF( tab_cntrl(27).EQ.1. ) ysinus =.true. 
+        IF( tab_cntrl(27).EQ.1. ) ysinus =.true.
       ENDIF
 
@@ -837 +837 @@
       real :: omega2(klon,klev),rho(klon,klev+1)
 !      real :: ug(klev),vg(klev),fcoriolis
-      real :: phis(klon) 
+      real :: phis(klon)
 
 !   Variables locales pour NetCDF:
@@ -893 +893 @@
 !       tab_cntrl(19) = preff
 !
-!    .....    parametres  pour le zoom      ......   
+!    .....    parametres  pour le zoom      ......
 
 !       tab_cntrl(20)  = clon
@@ -957 +957 @@
 !   passage d'un champ de la grille scalaire a la grille physique
 !=======================================================================
- 
+
 !-----------------------------------------------------------------------
 !   declarations:
 !   -------------
- 
+
       INTEGER im,jm,ngrid,nfield
       REAL pdyn(im,jm,nfield)
       REAL pfi(ngrid,nfield)
- 
+
       INTEGER i,j,ifield,ig
- 
+
 !-----------------------------------------------------------------------
 !   calcul:
 !   -------
- 
+
       DO ifield=1,nfield
 !   traitement des poles
@@ -978 +978 @@
             pdyn(i,jm,ifield)=pfi(ngrid,ifield)
          ENDDO
- 
+
 !   traitement des point normaux
          DO j=2,jm-1
@@ -986 +986 @@
          ENDDO
       ENDDO
- 
+
       RETURN
       END
- 
- 
+
+
 
       SUBROUTINE abort_gcm(modname, message, ierr)
- 
+
       USE IOIPSL
 !
@@ -1002 +1002 @@
 !         message = stuff to print
 !         ierr    = severity of situation ( = 0 normal )
- 
+
       character(len=*) modname
       integer ierr
       character(len=*) message
- 
+
       write(*,*) 'in abort_gcm'
       call histclo
@@ -1084 +1084 @@
       RETURN
       END
- 
+
       SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
       IMPLICIT NONE
@@ -1090 +1090 @@
 !   passage d'un champ de la grille scalaire a la grille physique
 !=======================================================================
- 
+
 !-----------------------------------------------------------------------
 !   declarations:
 !   -------------
- 
+
       INTEGER im,jm,ngrid,nfield
       REAL pdyn(im,jm,nfield)
       REAL pfi(ngrid,nfield)
- 
+
       INTEGER j,ifield,ig
- 
+
 !-----------------------------------------------------------------------
 !   calcul:
 !   -------
- 
+
       IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)                          &
      &    STOP 'probleme de dim'
@@ -1110 +1110 @@
       CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
       CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
- 
+
 !   traitement des point normaux
       DO ifield=1,nfield
@@ -1118 +1118 @@
          ENDDO
       ENDDO
- 
+
       RETURN
       END
- 
+
       SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
- 
+
 !    Ancienne version disvert dont on a modifie nom pour utiliser
 !    le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
@@ -1131 +1131 @@
 !
       USE dimensions_mod, ONLY: iim, jjm, llm, ndm
+USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
+          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
 IMPLICIT NONE
 
 
-      include "paramet.h"
+
 !
 !=======================================================================
@@ -1160 +1162 @@
       INTEGER np,ierr
       REAL pi,x
- 
+
 !-----------------------------------------------------------------------
 !
       pi=2.*ASIN(1.)
- 
+
       OPEN(99,file='sigma.def',status='old',form='formatted',                   &
      &   iostat=ierr)
- 
+
 !-----------------------------------------------------------------------
 !   cas 1 on lit les options dans sigma.def:
 !   ----------------------------------------
- 
+
       IF (ierr.eq.0) THEN
- 
+
       print*,'WARNING!!! on lit les options dans sigma.def'
       READ(99,*) deltaz
@@ -1184 +1186 @@
       alpha=deltaz/(llm*h)
 !
- 
+
        DO 1  l = 1, llm
        dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*                    &
@@ -1190 +1192 @@
      &            (1.-l/FLOAT(llm))*delta )
    1   CONTINUE
- 
+
        sig(1)=1.
        DO 101 l=1,llm-1
@@ -1196 +1198 @@
 101    CONTINUE
        sig(llm+1)=0.
- 
+
        DO 2  l = 1, llm
        dsig(l) = sig(l)-sig(l+1)
    2   CONTINUE
 !
- 
+
       ELSE
 !-----------------------------------------------------------------------
 !   cas 2 ancienne discretisation (LMD5...):
 !   ----------------------------------------
- 
+
       PRINT*,'WARNING!!! Ancienne discretisation verticale'
- 
+
       h=7.
       snorm  = 0.
@@ -1224 +1226 @@
          sig(l) = sig(l+1) + dsig(l)
       ENDDO
- 
+
       ENDIF
- 
- 
+
+
       DO l=1,llm
         nivsigs(l) = FLOAT(l)
       ENDDO
- 
+
       DO l=1,llmp1
         nivsig(l)= FLOAT(l)
       ENDDO
- 
+
 !
 !    ....  Calculs  de ap(l) et de bp(l)  ....
@@ -1243 +1245 @@
 !   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
 !
- 
+
       bp(llmp1) =   0.
- 
+
       DO l = 1, llm
 !c
 !cc    ap(l) = 0.
 !cc    bp(l) = sig(l)
- 
+
       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
       ap(l) = pa * ( sig(l) - bp(l) )
@@ -1256 +1258 @@
       ENDDO
       ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
- 
+
       PRINT *,' BP '
       PRINT *,  bp
       PRINT *,' AP '
       PRINT *,  ap
- 
+
       DO l = 1, llm
        dpres(l) = bp(l) - bp(l+1)
        presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
       ENDDO
- 
+
       PRINT *,' PRESNIVS '
       PRINT *,presnivs
- 
+
       RETURN
       END
@@ -1299 +1301 @@
 !   Schema amont pour l'advection verticale en 1D
 !   w est la vitesse verticale dp/dt en Pa/s
-!   Traitement en volumes finis 
+!   Traitement en volumes finis
 !   d / dt ( zm q ) = delta_z ( omega q )
 !   d / dt ( zm ) = delta_z ( omega )
@@ -1327 +1329 @@
       zwq(llm+1)=0.
       zw(llm+1)=0.
- 
+
       do l=1,llm
          qold=q(l)
@@ -1334 +1336 @@
       enddo
 
- 
+
       return
       end
@@ -1343 +1345 @@
        SUBROUTINE advect_va(llm,omega,d_t_va,d_q_va,d_u_va,d_v_va,              &
      &                q,temp,u,v,play)
-!itlmd 
+!itlmd
 !----------------------------------------------------------------------
-!   Calcul de l'advection verticale (ascendance et subsidence) de 
+!   Calcul de l'advection verticale (ascendance et subsidence) de
 !   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
-!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 
-!   sans WTG rajouter une advection horizontale 
-!----------------------------------------------------------------------  
+!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
+!   sans WTG rajouter une advection horizontale
+!----------------------------------------------------------------------
         implicit none
 INCLUDE "YOMCST.h"
@@ -1371 +1373 @@
      &       /(play(l)-play(l+1))
 
-        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))              
-
-        d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))              
-        d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))              
-
-        
+        d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))/(play(l)-play(l+1))
+
+        d_u_va(l)= -omgdown*(u(l)-u(l+1))/(play(l)-play(l+1))
+        d_v_va(l)= -omgdown*(v(l)-v(l+1))/(play(l)-play(l+1))
+
+
        elseif(l.eq.llm) then
         omgup=min(omega(l),0.0)
@@ -1387 +1389 @@
         d_u_va(l)= -omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
         d_v_va(l)= -omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
-       
+
        else
         omgup=min(omega(l),0.0)
@@ -1400 +1402 @@
         d_q_va(l,:)= -omgdown*(q(l,:)-q(l+1,:))                            &
      &              /(play(l)-play(l+1))-                                  &
-     &              omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l)) 
+     &              omgup*(q(l-1,:)-q(l,:))/(play(l-1)-play(l))
         d_u_va(l)= -omgdown*(u(l)-u(l+1))                                  &
      &              /(play(l)-play(l+1))-                                  &
-     &              omgup*(u(l-1)-u(l))/(play(l-1)-play(l)) 
+     &              omgup*(u(l-1)-u(l))/(play(l-1)-play(l))
         d_v_va(l)= -omgdown*(v(l)-v(l+1))                                  &
      &              /(play(l)-play(l+1))-                                  &
      &              omgup*(v(l-1)-v(l))/(play(l-1)-play(l))
-       
+
       endif
-         
+
       enddo
 !fin itlmd
@@ -1417 +1419 @@
        SUBROUTINE lstendH(llm,nqtot,omega,d_t_va,d_q_va,                        &
      &                q,temp,u,v,play)
-!itlmd 
+!itlmd
 !----------------------------------------------------------------------
-!   Calcul de l'advection verticale (ascendance et subsidence) de 
+!   Calcul de l'advection verticale (ascendance et subsidence) de
 !   temperature et d'humidite. Hypothese : ce qui rentre de l'exterieur
-!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou 
-!   sans WTG rajouter une advection horizontale 
-!----------------------------------------------------------------------  
+!   a les memes caracteristiques que l'air de la colonne 1D (WTG) ou
+!   sans WTG rajouter une advection horizontale
+!----------------------------------------------------------------------
         implicit none
 INCLUDE "YOMCST.h"
@@ -1648 +1650 @@
 !jyg<
 !   Formule pour q :
-!                         d_q = (1/tau) [rh_targ*qsat(T_new) - q] 
+!                         d_q = (1/tau) [rh_targ*qsat(T_new) - q]
 !
 !  Cette formule remplace d_q = (1/tau) [rh_targ - rh] qsat(T_new)
@@ -1743 +1745 @@
      &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas        &
      &         ,dth_mod_cas,hth_mod_cas,vth_mod_cas,mxcalc)
- 
+
        USE dimensions_mod, ONLY: iim, jjm, llm, ndm
 implicit none

LMDZ6/trunk/libf/phylmd/dyn1d/lmdz1d.F90

@@ -26 +26 @@
 
 
-INCLUDE "1DUTILS.h"
-INCLUDE "1Dconv.h"
+#include "1DUTILS.h"
+#include "1Dconv.h"
 
 !#endif

LMDZ6/trunk/libf/phylmd/dyn1d/paramet_mod_h.f90

@@ -1 @@
-link ../../dyn3d/paramet.h
+link ../../dyn3d/paramet_mod_h.f90