Changeset 4513 for LMDZ6/trunk/libf/phylmd/StratAer

Timestamp:

Apr 20, 2023, 9:57:22 AM (21 months ago)

Author:

tlurton

Message:

Strataer: a bug correction in strataer_mod.F90, and an implementation of injection duration for SAI case.

Location:

LMDZ6/trunk/libf/phylmd/StratAer

Files:

• strataer_mod.F90 (modified) (1 diff)
• traccoag_mod.F90 (modified) (5 diffs)

LMDZ6/trunk/libf/phylmd/StratAer/strataer_mod.F90

@@ -186 +186 @@
     !ENDIF
 
-    IF (flag_sulf_emit == 1) THEN
+    IF (flag_sulf_emit > 0) THEN
        CALL strataer_ponde_init
        WRITE(lunout,*) 'IN STRATAER INT2 END'

LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90

@@ -23 +23 @@
     USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     USE mod_phys_lmdz_para, only: gather, scatter
-    USE phys_cal_mod, ONLY : year_len, year_cur,mth_cur, day_cur, hour
+    USE phys_cal_mod, ONLY : year_len, mth_len, year_cur, mth_cur, day_cur, hour
     USE sulfate_aer_mod
     USE phys_local_var_mod, ONLY: stratomask
@@ -78 +78 @@
     REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
     REAL                                   :: dlat_loc
+    INTEGER                                :: injdur_sai          ! injection duration for SAI case [days]
+    INTEGER                                :: yr, is_bissext
 
     IF (is_mpi_root) THEN
@@ -229 +231 @@
     CASE(2) ! stratospheric aerosol injections (SAI)
 !
-      DO i=1,klon
-!       SAI standard scenario with continuous emission from 1 grid point at the equator
-!       SAI emission on single month
-!       SAI continuous emission o
+     ! Computing duration of SAI in days...
+     ! ... starting from 0...
+     injdur_sai = 0
+     ! ... then adding whole years from first to (n-1)th...
+     DO yr = year_emit_sai_start, year_emit_sai_end-1
+       ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
+       is_bissext = (MOD(yr,4)==0) .AND. (MOD(yr,100) /= 0 .OR. MOD(yr,400) == 0)
+       injdur_sai = injdur_sai+365+is_bissext
+     ENDDO
+     ! ... then subtracting part of the first year where no injection yet...
+     is_bissext = (MOD(year_emit_sai_start,4)==0) .AND. (MOD(year_emit_sai_start,100) /= 0 .OR. MOD(year_emit_sai_start,400) == 0)
+     SELECT CASE(mth_emit_sai_start)
+     CASE(2)
+        injdur_sai = injdur_sai-31
+     CASE(3)
+        injdur_sai = injdur_sai-31-28-is_bissext
+     CASE(4)
+        injdur_sai = injdur_sai-31-28-is_bissext-31
+     CASE(5)
+        injdur_sai = injdur_sai-31-28-is_bissext-31-30
+     CASE(6)
+        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31
+     CASE(7)
+        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30
+     CASE(8)
+        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31
+     CASE(9)
+        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31
+     CASE(10)
+        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30
+     CASE(11)
+        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31
+     CASE(12)
+        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31-30
+     END SELECT
+     injdur_sai = injdur_sai-day_emit_sai_start+1
+     ! ... then adding part of the n-th year 
+     is_bissext = (MOD(year_emit_sai_end,4)==0) .AND. (MOD(year_emit_sai_end,100) /= 0 .OR. MOD(year_emit_sai_end,400) == 0)
+     SELECT CASE(mth_emit_sai_end)
+     CASE(2)
+        injdur_sai = injdur_sai+31
+     CASE(3)
+        injdur_sai = injdur_sai+31+28+is_bissext
+     CASE(4)
+        injdur_sai = injdur_sai+31+28+is_bissext+31
+     CASE(5)
+        injdur_sai = injdur_sai+31+28+is_bissext+31+30
+     CASE(6)
+        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31
+     CASE(7)
+        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30
+     CASE(8)
+        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31
+     CASE(9)
+        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31
+     CASE(10)
+        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30
+     CASE(11)
+        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31
+     CASE(12)
+        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31+30
+     END SELECT
+     injdur_sai = injdur_sai+day_emit_sai_end
+     ! A security: are SAI dates of injection consistent?
+     IF (injdur_sai <= 0) THEN
+        CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.',1)
+     ENDIF
+     ! Injection in itself
+     IF (( year_emit_sai_start <= year_cur ) &
+        .AND. ( year_cur <= year_emit_sai_end ) &
+        .AND. ( mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur ) &
+        .AND. ( mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end ) &
+        .AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) &
+        .AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN
+
+     DO i=1,klon
         dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
         IF  ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. &
@@ -280 +354 @@
           DO k=1, klev
             ! stretch emission over whole year (360d)
-            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400.  
+            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(injdur_sai)/86400.  
             tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
             budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
@@ -295 +369 @@
         ENDIF ! emission grid cell
       ENDDO ! klon loop
+
+     ENDIF ! Condition over injection dates
 
     CASE(3) ! --- SAI injection over a single band of longitude and between