Changeset 1305 for trunk/LMDZ.VENUS/libf/phyvenus/physiq.F

Timestamp:

Jul 8, 2014, 2:21:53 PM (10 years ago)

Author:

slebonnois

Message:

SL: VENUS PHOTOCHEMISTRY. Needs Lapack (see arch files...)

File:

• trunk/LMDZ.VENUS/libf/phyvenus/physiq.F (modified) (8 diffs)

trunk/LMDZ.VENUS/libf/phyvenus/physiq.F

@@ -56 +56 @@
       USE ioipsl
 !      USE histcom ! not needed; histcom is included in ioipsl
+      USE chemparam_mod
       USE infotrac
       USE control_mod
@@ -65 +66 @@
       USE iophy
       use cpdet_mod, only: cpdet, t2tpot
+      use ieee_arithmetic
       IMPLICIT none
 c======================================================================
@@ -198 +200 @@
       REAL uq(klon) ! integr. verticale du transport zonal de l'eau
 c
+      REAL Fsedim(klon,klev+1)  ! Flux de sedimentation (kg.m-2)
 
 c======================================================================
@@ -285 +288 @@
       REAL u_seri(klon,klev), v_seri(klon,klev)
 c
-      REAL tr_seri(klon,klev,nqmax)
-      REAL d_tr(klon,klev,nqmax)
+      REAL :: tr_seri(klon,klev,nqmax)
+      REAL :: d_tr(klon,klev,nqmax)
+
+c Variables tendance sedimentation
+
+      REAL :: d_tr_sed(klon,klev,2)
+      REAL :: d_tr_ssed(klon)
 c
 c pour ioipsl
@@ -530 +538 @@
       ENDDO
 
+c---------
+c       Ecriture fichier initialisation
+c       PRINT*,'Ecriture Initial_State.csv'
+c       OPEN(88,file='Trac_Point.csv',
+c     & form='formatted')
+c---------
+     
+c---------
+c       Initialisation des parametres des nuages
+c===============================================
+     
+      if ((nlon .EQ. 1) .AND. ok_cloud) then
+        PRINT*,'Open profile_cloud_parameters.csv'
+        OPEN(66,file='profile_cloud_parameters.csv',
+     &   form='formatted')
+      endif
+
+      if ((nlon .EQ. 1) .AND. ok_sedim) then
+        PRINT*,'Open profile_cloud_sedim.csv'
+        OPEN(77,file='profile_cloud_sedim.csv',
+     &   form='formatted')
+      endif
+           
+      if ((nlon .GT. 1) .AND. ok_chem) then
+c !!! DONC 3D !!!
+        CALL chemparam_ini() 
+      endif
          
+      if ((nlon .GT. 1) .AND. ok_cloud) then
+c !!! DONC 3D !!!
+        CALL cloud_ini(nlon,nlev)
+      endif
+        
       ENDIF ! debut
 c====================================================================
@@ -704 +744 @@
        if (tr_scheme.eq.1) then
 ! Case 1: pseudo-chemistry with relaxation toward fixed profile
+
          call phytrac_relax (debut,lafin,nqmax,
      I                   nlon,nlev,dtime,pplay,
@@ -713 +754 @@
 ! However, the variable 'source' could be used in physiq
 ! so the call to phytrac_emiss could be to initialise it.
+
          call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
      I                   debut,lafin,nqmax,
@@ -718 +760 @@
      I                   rlatd,rlond,
      O                   tr_seri)
-       elseif (tr_scheme.eq.3) then
-! Case 3: Full chemistry
-!        call phytrac_chem ( ?? )
-         print*,"Chemistry not yet implemented..." 
-         print*,"See Aurelien Stolzenbach" 
-       endif
-      endif
+
+       elseif (tr_scheme.eq.3) then  ! identical to ok_chem.or.ok_cloud
+! Case 3: Full chemistry and/or clouds
+
+         call phytrac_chimie(                  
+     I             debut,
+     I             gmtime,
+     I             nqmax,
+     I             nlon,
+     I             rlatd,
+     I             rlond,
+     I             nlev,
+     I             dtime,
+     I             t_seri,pplay,
+     O             tr_seri,
+     O             NBRTOT,
+     O             WH2SO4,
+     O             rho_droplet)
+
+c        CALL WriteField_phy('Pression',pplay,nlev)
+c        CALL WriteField_phy('PressionBnd',paprs,nlev+1)
+c        CALL WriteField_phy('Temp',t_seri,nlev)
+c        IF (ok_cloud) THEN
+c          CALL WriteField_phy('NBRTOT',NBRTOT,nlev)
+c        ENDIF
+c        CALL WriteField_phy('SAl',tr_seri(:,:,i_h2so4liq),nlev)
+c        CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev)
+
+         if (ok_sedim) then 
+
+           CALL new_cloud_sedim(
+     I               klon,
+     I               nlev,
+     I               dtime,
+     I               pplay,
+     I               paprs,
+     I               t_seri, 
+     I               WH2SO4,
+     I               tr_seri,
+     I               nqmax,
+     I               NBRTOT,
+     I               rho_droplet,
+     O               Fsedim,
+     O               d_tr_sed,
+     O               d_tr_ssed)
+
+          DO k = 1, klev
+           DO i = 1, klon
+     
+c        WRITE(88,"(11(e15.8,','))") pplay(5,25),
+c     &  t_seri(5,25),tr_seri(5,25,i_h2oliq),
+c     &  tr_seri(5,25,i_h2o),tr_seri(5,25,i_h2so4liq),
+c     &  tr_seri(5,25,i_h2so4),NBRTOT(5,25),WH2SO4(5,25),
+c     &  Fsedim(5,25),d_tr_sed(5,25,1),d_tr_sed(5,25,2)
+
+c--------------------
+c   Ce test est necessaire pour eviter Xliq=NaN   
+        IF (ieee_is_nan(d_tr_sed(i,k,1)).OR.
+     &  ieee_is_nan(d_tr_sed(i,k,2))) THEN
+        PRINT*,'sedim NaN PROBLEM'
+        PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k)
+        PRINT*,'lat-lon',i,'level',k,'dtime',dtime
+        PRINT*,'NBRTOT',NBRTOT(i,k),'F_sed',Fsedim(i,k)
+        PRINT*,'==============================================='
+                d_tr_sed(i,k,:)=0.
+        ENDIF
+c--------------------
+
+        tr_seri(i,k,i_h2so4liq) = tr_seri(i,k,i_h2so4liq)+
+     &                            d_tr_sed(i,k,1)
+        tr_seri(i,k,i_h2oliq)   = tr_seri(i,k,i_h2oliq)+
+     &                            d_tr_sed(i,k,2)
+        d_tr_sed(i,k,:) = d_tr_sed(i,k,:) / dtime
+        Fsedim(i,k)     = Fsedim(i,k) / dtime
+     
+           ENDDO
+          ENDDO
+      
+        Fsedim(:,klev+1) = 0.
+
+         endif ! ok_sedim
+
+       endif   ! tr_scheme
+      endif    ! iflag_trac
 
 c