Changeset 4171 for LMDZ6/branches/LMDZ-ECRAD/libf/dyn3dmem/iniacademic_loc.F90

Timestamp:

Jun 17, 2022, 4:24:49 PM (2 years ago)

Author:

lguez

Message:

Sync latest trunk changes to branch LMDZ-ECRAD.

Location:

LMDZ6/branches/LMDZ-ECRAD

Files:

• . (modified) (1 prop)
• libf/dyn3dmem/iniacademic_loc.F90 (modified) (11 diffs)

LMDZ6/branches/LMDZ-ECRAD/libf/dyn3dmem/iniacademic_loc.F90

@@ -5 +5 @@
 
   USE filtreg_mod, ONLY: inifilr
+  USE infotrac,    ONLY: nqtot, niso, tnat, alpha_ideal, iqIsoPha, tracers
+  USE control_mod, ONLY: day_step,planet_type
   use exner_hyb_m, only: exner_hyb
   use exner_milieu_m, only: exner_milieu
-  USE infotrac, ONLY: nqtot,niso_possibles,ok_isotopes,iqpere,ok_iso_verif,tnat,alpha_ideal, &
-        & iqiso,phase_num,iso_indnum,iso_num,zone_num
-  USE control_mod, ONLY: day_step,planet_type
   USE parallel_lmdz, ONLY: ijb_u, ije_u, ijb_v, ije_v
 #ifdef CPP_IOIPSL
@@ -23 +22 @@
   USE temps_mod, ONLY: annee_ref, day_ini, day_ref
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
-
+  USE readTracFiles_mod, ONLY: addPhase
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -67 +66 @@
   real tetastrat ! potential temperature in the stratosphere, in K
   real tetajl(jjp1,llm)
-  INTEGER i,j,l,lsup,ij
+  INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
 
   REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
@@ -73 +72 @@
   LOGICAL ok_geost             ! Initialisation vent geost. ou nul
   LOGICAL ok_pv                ! Polar Vortex
-  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
+  REAL phi_pv,dphi_pv,gam_pv,tetanoise   ! Constantes pour polar vortex 
 
   real zz,ran1
@@ -112 +111 @@
   ztot0      = 0.
   stot0      = 0.
-  ang0       = 0.      
+  ang0       = 0.
 
   if (llm == 1) then
@@ -122 +121 @@
   CALL inigeom
   CALL inifilr
+
+  ! Initialize pressure and mass field if read_start=.false.
+  IF (.NOT. read_start) THEN
+    ! allocate global fields:
+!    allocate(vcov_glo(ip1jm,llm))
+    allocate(ucov_glo(ip1jmp1,llm))
+    allocate(teta_glo(ip1jmp1,llm))
+    allocate(ps_glo(ip1jmp1))
+    allocate(masse_glo(ip1jmp1,llm))
+    allocate(phis_glo(ip1jmp1))
+
+     ! surface pressure
+     if (iflag_phys>2) then
+        ! specific value for CMIP5 aqua/terra planets
+        ! "Specify the initial dry mass to be equivalent to
+        !  a global mean surface pressure (101325 minus 245) Pa."
+        ps_glo(:)=101080.  
+     else
+        ! use reference surface pressure
+        ps_glo(:)=preff
+     endif
+
+     ! ground geopotential
+     phis_glo(:)=0.
+     CALL pression ( ip1jmp1, ap, bp, ps_glo, p       )
+     if (pressure_exner) then
+       CALL exner_hyb( ip1jmp1, ps_glo, p, pks, pk )
+     else
+       call exner_milieu(ip1jmp1,ps_glo,p,pks,pk)
+     endif
+     CALL massdair(p,masse_glo)
+  ENDIF
 
   if (llm == 1) then
@@ -172 +203 @@
      gam_pv=4.              ! -dT/dz vortex (in K/km)
      CALL getin('gam_pv',gam_pv)
+     tetanoise=0.005
+     CALL getin('tetanoise',tetanoise)
 
      ! 2. Initialize fields towards which to relax
@@ -224 +257 @@
      ! 3. Initialize fields (if necessary)
      IF (.NOT. read_start) THEN
-       ! allocate global fields:
-!       allocate(vcov_glo(ip1jm,llm))
-       allocate(ucov_glo(ip1jmp1,llm))
-       allocate(teta_glo(ip1jmp1,llm))
-       allocate(ps_glo(ip1jmp1))
-       allocate(masse_glo(ip1jmp1,llm))
-       allocate(phis_glo(ip1jmp1))
-
-        ! surface pressure
-        if (iflag_phys>2) then
-           ! specific value for CMIP5 aqua/terra planets
-           ! "Specify the initial dry mass to be equivalent to
-           !  a global mean surface pressure (101325 minus 245) Pa."
-           ps_glo(:)=101080.  
-        else
-           ! use reference surface pressure
-           ps_glo(:)=preff
-        endif
-        
-        ! ground geopotential
-        phis_glo(:)=0.
-
-        CALL pression ( ip1jmp1, ap, bp, ps_glo, p       )
-        if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps_glo, p, pks, pk )
-        else
-          call exner_milieu(ip1jmp1,ps_glo,p,pks,pk)
-        endif
-        CALL massdair(p,masse_glo)
-
         ! bulk initialization of temperature
-        teta_glo(:,:)=tetarappel(:,:)
-
+        IF (iflag_phys>10000) THEN
+        ! Particular case to impose a constant temperature T0=0.01*iflag_phys
+           teta_glo(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp)
+        ELSE
+           teta_glo(:,:)=tetarappel(:,:)
+        ENDIF
         ! geopotential
         CALL geopot(ip1jmp1,teta_glo,pk,pks,phis_glo,phi)
@@ -271 +278 @@
         if (planet_type=="earth") then
            ! Earth: first two tracers will be water
-
-           do i=1,nqtot
-              if (i == 1) q(ijb_u:ije_u,:,i)=1.e-10
-              if (i == 2) q(ijb_u:ije_u,:,i)=1.e-15
-              if (i.gt.2) q(ijb_u:ije_u,:,i)=0.
+           do iq=1,nqtot
+              q(ijb_u:ije_u,:,iq)=0.
+              IF(tracers(iq)%name == addPhase('H2O', 'g')) q(ijb_u:ije_u,:,iq)=1.e-10
+              IF(tracers(iq)%name == addPhase('H2O', 'l')) q(ijb_u:ije_u,:,iq)=1.e-15
 
               ! CRisi: init des isotopes
               ! distill de Rayleigh très simplifiée
-              if (ok_isotopes) then
-                if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then          
-                   q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqpere(i))       &
-      &                  *tnat(iso_num(i))                             &
-      &                  *(q(ijb_u:ije_u,:,iqpere(i))/30.e-3)                              &
-     &                   **(alpha_ideal(iso_num(i))-1)
-                endif                
-                if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then
-                  q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqiso(iso_indnum(i),phase_num(i)))
-                endif
-              endif !if (ok_isotopes) then
-
+              iName    = tracers(iq)%iso_iName
+              if (niso <= 0 .OR. iName <= 0) CYCLE
+              iPhase   = tracers(iq)%iso_iPhase
+              iqParent = tracers(iq)%iqParent
+              IF(tracers(iq)%iso_iZone == 0) THEN
+                 q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat(iName) &
+                                     *(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+              ELSE
+                 q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
+              END IF
            enddo
         else
@@ -296 +300 @@
         endif ! of if (planet_type=="earth")
 
-        if (ok_iso_verif) then
-           call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc')
-        endif !if (ok_iso_verif) then
+        call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc')
 
         ! add random perturbation to temperature
@@ -306 +308 @@
         do l=1,llm
            do ij=iip2,ip1jm
-              teta_glo(ij,l)=teta_glo(ij,l)*(1.+0.005*ran1(idum))
+              teta_glo(ij,l)=teta_glo(ij,l)*(1.+tetanoise*ran1(idum))
            enddo
         enddo