Changeset 5281 for LMDZ6/trunk/libf/dyn3d/guide_mod.f90

Timestamp:

Oct 28, 2024, 11:17:48 AM (3 weeks ago)

Author:

abarral

Message:

Turn comgeom.h comgeom2.h into modules

File:

• LMDZ6/trunk/libf/dyn3d/guide_mod.f90 (modified) (43 diffs)

LMDZ6/trunk/libf/dyn3d/guide_mod.f90

@@ -8 +8 @@
 !            F. Codron 01/09
 !=======================================================================
-
-  USE getparam, only: ini_getparam, fin_getparam, getpar
+    USE getparam, only: ini_getparam, fin_getparam, getpar
   USE Write_Field
   use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, &
@@ -92 +91 @@
 
 ! ---------------------------------------------
-! Lecture des parametres:  
+! Lecture des parametres:
 ! ---------------------------------------------
     call ini_getparam("nudging_parameters_out.txt")
@@ -172 +171 @@
 
     call fin_getparam
-    
+
 ! ---------------------------------------------
 ! Determination du nombre de niveaux verticaux
@@ -230 +229 @@
 
 
-    endif 
+    endif
     error=nf90_inq_dimid(ncidpl,'LEVEL',rid)
     IF (error.NE.NF90_NOERR) error=nf90_inq_dimid(ncidpl,'PRESSURE',rid)
@@ -237 +236 @@
     ENDIF
     error=nf90_inquire_dimension(ncidpl,rid,len=nlevnc)
-    write(*,*)trim(modname)//' : number of vertical levels nlevnc', nlevnc 
+    write(*,*)trim(modname)//' : number of vertical levels nlevnc', nlevnc
     rcod = nf90_close(ncidpl)
 
@@ -264 +263 @@
     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
     alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
-    
+
     IF (guide_u) THEN
         ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
@@ -288 +287 @@
         tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
     ENDIF
-     
+
     IF (guide_Q) THEN
         ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
@@ -361 +360 @@
     USE comconst_mod, ONLY: cpp, dtvr, daysec,kappa
     USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner
- 
+
     USE dimensions_mod, ONLY: iim, jjm, llm, ndm
 USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
@@ -391 +390 @@
     REAL          :: tau,reste ! position entre 2 etats de guidage
     REAL, SAVE    :: factt ! pas de temps en fraction de jour
-    
+
     INTEGER       :: l
     CHARACTER(LEN=20) :: modname="guide_main"
@@ -402 +401 @@
     IF (first) THEN
         first=.FALSE.
-        CALL guide_init 
+        CALL guide_init
         itau_test=1001
         step_rea=1
         count_no_rea=0
 ! Calcul des constantes de rappel
-        factt=dtvr*iperiod/daysec 
+        factt=dtvr*iperiod/daysec
         call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
         call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
@@ -421 +420 @@
             enddo
         endif
-! ini_anal: etat initial egal au guidage        
+! ini_anal: etat initial egal au guidage
         IF (ini_anal) THEN
             CALL guide_interp(ps,teta)
@@ -464 +463 @@
             itau
               abort_message='stopped'
-              CALL abort_gcm(modname,abort_message,1) 
+              CALL abort_gcm(modname,abort_message,1)
           ELSE
               IF (guide_v) vnat1=vnat2
@@ -503 +502 @@
 
 !-----------------------------------------------------------------------
-!   Ajout des champs de guidage 
+!   Ajout des champs de guidage
 !-----------------------------------------------------------------------
 ! Sauvegarde du guidage?
-    f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)  
+    f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
     IF (f_out) THEN
       ! compute pressures at layer interfaces
@@ -522 +521 @@
       CALL guide_out("SP",jjp1,llm,p(:,1:llm))
     ENDIF
-    
+
     if (guide_u) then
         if (guide_add) then
@@ -528 +527 @@
         else
            f_add=(1.-tau)*ugui1+tau*ugui2-ucov
-        endif 
+        endif
         if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
@@ -542 +541 @@
         else
            f_add=(1.-tau)*tgui1+tau*tgui2-teta
-        endif 
+        endif
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
@@ -554 +553 @@
         else
            f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2-ps
-        endif 
+        endif
         if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
         CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
@@ -568 +567 @@
         else
            f_add=(1.-tau)*qgui1+tau*qgui2-q
-        endif 
+        endif
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
@@ -580 +579 @@
         else
            f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2-vcov
-        endif 
+        endif
         if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
         CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
@@ -599 +598 @@
     INTEGER,                      INTENT(IN)    :: hsize
     INTEGER,                      INTENT(IN)    :: vsize
-    REAL, DIMENSION(hsize),       INTENT(IN)    :: alpha 
+    REAL, DIMENSION(hsize),       INTENT(IN)    :: alpha
     REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
 
@@ -615 +614 @@
 
     USE comconst_mod, ONLY: pi
-    
+
     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
+    USE comgeom_mod_h
 IMPLICIT NONE
 
 
 
-    INCLUDE "comgeom.h"
-    
+
     ! input/output variables
     INTEGER,                           INTENT(IN)    :: typ
@@ -665 +664 @@
                 fieldm(j,l)=fieldm(j,l)+field(ij,l)
             ENDDO
-        ENDDO 
+        ENDDO
         fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
     ! Compute forcing
@@ -680 +679 @@
 !=======================================================================
   SUBROUTINE guide_interp(psi,teta)
-  
+
   use exner_hyb_m, only: exner_hyb
   use exner_milieu_m, only: exner_milieu
@@ -686 +685 @@
   use comvert_mod, only: preff, pressure_exner, bp, ap
   USE dimensions_mod, ONLY: iim, jjm, llm, ndm
-USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
+  USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
           ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
+  USE comgeom2_mod_h
 IMPLICIT NONE
-
-
-
-  include "comgeom2.h"
-
   REAL, DIMENSION (iip1,jjp1),     INTENT(IN) :: psi ! Psol gcm
   REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
@@ -702 +697 @@
   REAL, DIMENSION (iip1,jjp1,llm)    :: plunc,plsnc !niveaux pression modele
   REAL, DIMENSION (iip1,jjm,llm)     :: plvnc       !niveaux pression modele
-  REAL, DIMENSION (iip1,jjp1,llmp1)  :: p           ! pression intercouches 
+  REAL, DIMENSION (iip1,jjp1,llmp1)  :: p           ! pression intercouches
   REAL, DIMENSION (iip1,jjp1,llm)    :: pls, pext   ! var intermediaire
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pbarx 
-  REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
+  REAL, DIMENSION (iip1,jjp1,llm)    :: pbarx
+  REAL, DIMENSION (iip1,jjm,llm)     :: pbary
   ! Variables pour fonction Exner (P milieu couche)
   REAL, DIMENSION (iip1,jjp1,llm)    :: pk
-  REAL, DIMENSION (iip1,jjp1)        :: pks    
+  REAL, DIMENSION (iip1,jjp1)        :: pks
   REAL                               :: prefkap,unskap
   ! Pression de vapeur saturante
   REAL, DIMENSION (ip1jmp1,llm)      :: qsat
   !Variables intermediaires interpolation
-  REAL, DIMENSION (iip1,jjp1,llm)    :: zu1,zu2 
+  REAL, DIMENSION (iip1,jjp1,llm)    :: zu1,zu2
   REAL, DIMENSION (iip1,jjm,llm)     :: zv1,zv2
-  
+
   INTEGER                            :: i,j,l,ij
   CHARACTER(LEN=20),PARAMETER :: modname="guide_interp"
-  
+
     write(*,*)trim(modname)//': interpolate nudging variables'
 ! -----------------------------------------------------------------
@@ -767 +762 @@
         endif
     endif
-    
+
 ! -----------------------------------------------------------------
-! Calcul niveaux pression modele 
+! Calcul niveaux pression modele
 ! -----------------------------------------------------------------
     CALL pression( ip1jmp1, ap, bp, psi, p )
@@ -839 +834 @@
                     ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
                 enddo
-                ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)    
-                ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)    
+                ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)
+                ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)
             enddo
             do i=1,iip1
@@ -850 +845 @@
         enddo
     ENDIF
-    
+
     IF (guide_T) THEN
         CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
@@ -869 +864 @@
                     enddo
                 ENDIF
-                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)    
-                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)    
+                tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)
+                tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)
             enddo
             do i=1,iip1
                 tgui1(i,l)=tgui1(1,l)
-                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l) 
+                tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
                 tgui2(i,l)=tgui2(1,l)
-                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l) 
+                tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
             enddo
         enddo
@@ -893 +888 @@
                     vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
                 enddo
-                vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)    
-                vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)    
+                vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)
+                vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)
             enddo
         enddo
     ENDIF
-    
+
     IF (guide_Q) THEN
         ! On suppose qu'on a la bonne variable dans le fichier de guidage:
@@ -911 +906 @@
                     qgui2(ij,l)=zu2(i,j,l)
                 enddo
-                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)    
-                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)    
+                qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)
+                qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)
             enddo
             do i=1,iip1
                 qgui1(i,l)=qgui1(1,l)
-                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l) 
+                qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
                 qgui2(i,l)=qgui2(1,l)
-                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l) 
+                qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
             enddo
         enddo
@@ -924 +919 @@
             CALL q_sat(iip1*jjp1*llm,teta*pk/cpp,plsnc,qsat)
             qgui1=qgui1*qsat*0.01 !hum. rel. en %
-            qgui2=qgui2*qsat*0.01 
+            qgui2=qgui2*qsat*0.01
         ENDIF
     ENDIF
@@ -937 +932 @@
     use comconst_mod, only: pi
     use serre_mod, only: clon, clat, grossismx, grossismy
-    
+
     USE dimensions_mod, ONLY: iim, jjm, llm, ndm
-USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
+    USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
           ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
-implicit none
-
-
-
-    include "comgeom2.h"
-
+    USE comgeom2_mod_h
+    implicit none
 ! input arguments :
     INTEGER, INTENT(IN) :: typ    ! u(2),v(3), ou scalaire(1)
@@ -953 +944 @@
     REAL, INTENT(IN)    :: taumin,taumax
 ! output arguments:
-    REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha 
-  
+    REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
+
 !  local variables:
     LOGICAL, SAVE               :: first=.TRUE.
@@ -1040 +1031 @@
         ENDIF
 ! Premier appel: calcul des aires min et max et de gamma.
-        IF (first) THEN 
+        IF (first) THEN
             first=.FALSE.
             ! coordonnees du centre du zoom
-            CALL coordij(clon,clat,ilon,ilat) 
+            CALL coordij(clon,clat,ilon,ilat)
             ! aire de la maille au centre du zoom
             dxdy_min=dxdys(ilon,ilat)
@@ -1067 +1058 @@
               endif
               gamma=log(0.5)/log(gamma)
-              if (gamma4) then 
+              if (gamma4) then
                 gamma=min(gamma,4.)
               endif
@@ -1168 +1159 @@
              ENDIF
              rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
-             IF (rcode.NE.NF90_NOERR) THEN 
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no PRES variable in file P.nc'
               CALL abort_gcm(modname,abort_message,1)
@@ -1221 +1212 @@
              write(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv
              if (ncidpl.eq.-99) ncidpl=ncidv
-             
+
              status=NF90_INQ_DIMID(ncidv, "LONV", dimid)
              status=NF90_INQUIRE_DIMENSION(ncidv,dimid,namedim,lendim)
-             
+
                 IF (lendim .NE. iip1) THEN
                 abort_message='dimension LONV different from iip1 in v.nc'
@@ -1237 +1228 @@
                 CALL abort_gcm(modname,abort_message,1)
              ENDIF
-        
+
          endif
 
@@ -1350 +1341 @@
      count(4)=1
 
-! Pression 
+! Pression
      if (guide_plevs.EQ.2) then
          status=NF90_GET_VAR(ncidp,varidp,pnat2,start,count)
@@ -1382 +1373 @@
            CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
          ENDIF
-         
+
      endif
 
@@ -1672 +1663 @@
 
   END SUBROUTINE guide_read2D
-  
+
 !=======================================================================
   SUBROUTINE guide_out(varname,hsize,vsize,field)
@@ -1680 +1671 @@
     use netcdf95, only: nf95_def_var, nf95_put_var
     use netcdf, only: nf90_float, nf90_def_var, nf90_put_var
-    
+
     USE dimensions_mod, ONLY: iim, jjm, llm, ndm
-USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
+    USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
           ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
+    USE comgeom2_mod_h
 IMPLICIT NONE
 
 
 
-    INCLUDE "comgeom2.h"
     
     ! Variables entree