Changeset 2359 for LMDZ5/trunk/libf/phylmd/physiq.F90

Timestamp:

Sep 4, 2015, 6:48:31 PM (9 years ago)

Author:

oboucher

Message:

Recoded diurnal cycle in physiq.F90
iflag_diurnal_cycle is the new flag
=0 for no diurnal cycle
=1 for old diurnal cycle
=2 for new diurnal cycle with updates
in between two calls to radiation according
to how rmu0 has changed over the previous
timestep of physiq

File:

• LMDZ5/trunk/libf/phylmd/physiq.F90 (modified) (11 diffs)

LMDZ5/trunk/libf/phylmd/physiq.F90

@@ -143 +143 @@
   !======================================================================
   ! Clef controlant l'activation du cycle diurne:
-  !cc      LOGICAL cycle_diurne
-  !cc      PARAMETER (cycle_diurne=.FALSE.)
+  ! en attente du codage des cles par Fred
+        INTEGER iflag_cycle_diurne
+        PARAMETER (iflag_cycle_diurne=1)
   !======================================================================
   ! Modele thermique du sol, a activer pour le cycle diurne:
@@ -589 +590 @@
   !
   REAL dist, rmu0(klon), fract(klon)
-  REAL zdtime, zlongi
+  REAL zrmu0(klon), zfract(klon), swradcorr(klon)
+  REAL zdtime, zdtime1, zdtime2, zlongi
   !
   REAL qcheck
@@ -1070 +1072 @@
      print*,'iflag_coupl,iflag_clos,iflag_wake', &
           iflag_coupl,iflag_clos,iflag_wake
-     print*,'CYCLE_DIURNE', cycle_diurne
+     print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
      !
      IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
@@ -1122 +1124 @@
      ENDIF
 
-     !IM cf. AM 081204 BEG
-     PRINT*,'cycle_diurne3 =',cycle_diurne
-     !IM cf. AM 081204 END
-     !
      CALL printflag( tabcntr0,radpas,ok_journe, &
           ok_instan, ok_region )
@@ -1149 +1147 @@
      ENDIF
      !
-     IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN 
+     IF (dtime*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN 
         WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
         WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
@@ -1821 +1819 @@
   ! non nul aux poles.
   IF (abs(solarlong0-1000.)<1.e-4) then
-     call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
+     call zenang_an(iflag_cycle_diurne.GE.1,jH_cur,rlat,rlon,rmu0,fract)
   ELSE
-     !  Avec ou sans cycle diurne
-     IF (cycle_diurne) THEN
+  ! recode par Olivier Boucher en sept 2015
+     SELECT CASE (iflag_cycle_diurne)
+     CASE(0)  
+     !  Sans cycle diurne
+        CALL angle(zlongi, rlat, fract, rmu0)
+     CASE(1)  
+     !  Avec cycle diurne sans application des poids
+     !  bit comparable a l ancienne formulation cycle_diurne=true
+     !  on integre entre gmtime et gmtime+radpas
         zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
-        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
-     ELSE
-        CALL angle(zlongi, rlat, fract, rmu0)
-     ENDIF
+        CALL zenang(zlongi,jH_cur,0.0,zdtime,rlat,rlon,rmu0,fract)
+     CASE(2)  
+     !  Avec cycle diurne sans application des poids
+     !  On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
+     !  Comme cette routine est appele a tous les pas de temps de la physique
+     !  meme si le rayonnement n'est pas appele je remonte en arriere les
+     !  radpas-1 pas de temps suivant. Petite ruse pour prendre en compte le
+     !  premier pas de temps la physique ou itaprad=0 
+        zdtime1=dtime*REAL(-MOD(itaprad,4)-1)      
+        zdtime2=dtime*REAL(radpas-MOD(itaprad,4)-1) 
+        CALL zenang(zlongi,jH_cur,zdtime1,zdtime2,rlat,rlon,rmu0,fract)
+     END SELECT 
   ENDIF
-
-  ! AI Janv 2014 
-  do i = 1, klon
-     if (fract(i).le.0.) then
-        JrNt(i)=0.
-     else
-        JrNt(i)=1.
-     endif
-  enddo
 
   if (mydebug) then
@@ -3562 +3566 @@
              topsw0,toplw0,solsw0,sollw0, &
              lwdn0, lwdn, lwup0, lwup,  &
-             swdn0, swdn, swup0, swup, &
+             zswdn0, zswdn, zswup0, zswup, &
              ok_ade, ok_aie, &
              tau_aero, piz_aero, cg_aero, &
@@ -3608 +3612 @@
              topsw0,toplw0,solsw0,sollw0, &
              lwdn0, lwdn, lwup0, lwup,  &
-             swdn0, swdn, swup0, swup, &
+             zswdn0, zswdn, zswup0, zswup, &
              topswad_aero, solswad_aero, &
              topswai_aero, solswai_aero, &
@@ -3664 +3668 @@
                    topsw0p,toplw0p,solsw0p,sollw0p, &
                    lwdn0p, lwdnp, lwup0p, lwupp,  &
-                   swdn0p, swdnp, swup0p, swupp, &
+                   zswdn0p, zswdnp, zswup0p, zswupp, &
                    topswad_aerop, solswad_aerop, &
                    topswai_aerop, solswai_aerop, &
@@ -3706 +3710 @@
      swup=0.    ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
      swup0=0.
-     swdn=0.
-     swdn0=0.
+     zswdn=0.
+     zswdn0=0.
      lwup=0.
      lwup0=0.
@@ -3715 +3719 @@
 
   !
+  ! Calculer les poids a appliquer sur le SW 
+  ! sortie JrNt = jour-nuit 
+  ! recode par Olivier Boucher en sept 2015
+  ! 
+
+  SELECT CASE (iflag_cycle_diurne)
+  ! 
+  CASE(0)
+  !  Sans cycle diurne
+     swradcorr=1.0
+     JrNt = 1.0
+  CASE(1)
+  !  Avec cycle diurne sans les poids
+     swradcorr=1.0
+     JrNt=0.0
+     WHERE (fract.GT.0.0) JrNt=1.0
+  CASE(2)
+  !  Avec cycle diurne et les poids
+     zdtime1=-dtime !--on corrige le rayonnement pour representer le
+     zdtime2=0.0    !--pas de temps de la physique qui se termine
+     CALL zenang(zlongi,jH_cur,zdtime1,zdtime2,rlat,rlon,zrmu0,zfract)
+     swradcorr=0.0
+     WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) swradcorr=zfract/fract*zrmu0/rmu0
+     JrNt=0.0
+     WHERE (zfract.GT.0.0) JrNt=1.0 
+  END SELECT
+
+  !
+  ! Corriger les flux SW pour le cycle diurne ameliore
+  ! recode par Olivier Boucher en sept 2015
+  !
+
+  DO k=1, klev+1
+    swdn0(:,k)=swradcorr(:)*zswdn0(:,k)
+    swdn(:,k) =swradcorr(:)*zswdn(:,k)
+    swup0(:,k)=swradcorr(:)*zswup0(:,k)
+    swup(:,k) =swradcorr(:)*zswup(:,k)
+  ENDDO
+  if (ok_4xCO2atm) then
+  DO k=1, klev+1
+    swdn0p(:,k)=swradcorr(:)*zswdn0p(:,k)
+    swdnp(:,k) =swradcorr(:)*zswdnp(:,k)
+    swup0p(:,k)=swradcorr(:)*zswup0p(:,k)
+    swupp(:,k) =swradcorr(:)*zswupp(:,k)
+  ENDDO
+  endif
+
+  !
   ! Ajouter la tendance des rayonnements (tous les pas)
   !
-  d_t_swr(:,:)=heat(:,:)*dtime/RDAY
-  d_t_lwr(:,:)=-cool(:,:)*dtime/RDAY
-  d_t_sw0(:,:)=heat0(:,:)*dtime/RDAY
-  d_t_lw0(:,:)=-cool0(:,:)*dtime/RDAY
+  
+  DO k=1, klev
+    d_t_swr(:,k)=swradcorr(:)*heat(:,k)*dtime/RDAY
+    d_t_sw0(:,k)=swradcorr(:)*heat0(:,k)*dtime/RDAY
+    d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
+    d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
+  ENDDO
+
   CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy)
   CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy)