Changeset 1975 for branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90

Timestamp:

2010-06-28T19:22:14+02:00 (14 years ago)

Author:

mlelod

Message:

ticket: #663 MLF: first part

File:

• branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90 (modified) (10 diffs)

branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/traqsr.F90

@@ -47 +47 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
    
-   INTEGER                   ::   nksr    ! levels below which the light cannot penetrate (depth larger than 391 m)
+   INTEGER , PUBLIC          ::   nksr    ! levels below which the light cannot penetrate (depth larger than 391 m)
    REAL(wp), DIMENSION(3,61) ::   rkrgb   ! tabulated attenuation coefficients for RGB absorption
 
@@ -97 +97 @@
       REAL(wp) ::   zchl, zcoef, zsi0r   ! temporary scalars
       REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
-      REAL(wp) ::   zqsr                 !    -         -
       REAL(wp), DIMENSION(jpi,jpj)     ::   zekb, zekg, zekr            ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze0, ze1 , ze2, ze3, zea    ! 3D workspace
@@ -123 +122 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
 !!gm  how to stecify the mean of time step here : TOP versus OPA time stepping strategy not obvious
-                  ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
+                  qsr_trd_hc_n(ji,jj,jk) = ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
                END DO
             END DO
@@ -155 +154 @@
                zsi0r = 1.e0 / rn_si0
                zcoef = ( 1. - rn_abs ) / 3.e0                           ! equi-partition in R-G-B
-!!gm bug !!! zqsr only a constant not an array
-               zqsr  = 0.5 * ( qsr_b(:,:) + qsr(:,:) )                  ! mean over 2 time steps
-               ze0(:,:,1) = rn_abs  * zqsr
-               ze1(:,:,1) = zcoef   * zqsr
-               ze2(:,:,1) = zcoef   * zqsr
-               ze3(:,:,1) = zcoef   * zqsr
-               zea(:,:,1) =           zqsr
+
+               ze0(:,:,1) = rn_abs  * qsr(:,:)
+               ze1(:,:,1) = zcoef   * qsr(:,:)
+               ze2(:,:,1) = zcoef   * qsr(:,:)
+               ze3(:,:,1) = zcoef   * qsr(:,:)
+               zea(:,:,1) =           qsr(:,:)
                !
                DO jk = 2, nksr+1                                     ! deeper values
@@ -183 +181 @@
                !
                DO jk = 1, nksr                                       ! compute and add qsr trend to ta
-                  ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+                  qsr_trd_hc_n(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) )
                END DO
                zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero
@@ -190 +188 @@
             ELSE                                                 !*  Constant Chlorophyll
                DO jk = 1, nksr
-                  ta(:,:,jk) = ta(:,:,jk) + etot3(:,:,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) )
+                  qsr_trd_hc_n(:,:,jk) = etot3(:,:,jk) * qsr(:,:)
                END DO
             ENDIF
@@ -201 +199 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ta(ji,jj,jk) = ta(ji,jj,jk) + etot3(ji,jj,jk) * 0.5 * ( qsr_b(:,:) + qsr(:,:) )
+                     qsr_trd_hc_n(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj)
                   END DO
                END DO
@@ -208 +206 @@
          ENDIF
          !
+      ENDIF
+
+      ! Add qsr trend to ta in all cases
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+         DO jk = 1, nksr
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ta(ji,jj,jk) = ta(ji,jj,jk) + qsr_trd_hc_n(ji,jj,jk) / fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE
+         DO jk = 1, nksr
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ta(ji,jj,jk) = ta(ji,jj,jk) + 0.5 * ( qsr_trd_hc_b(ji,jj,jk) + qsr_trd_hc_n(ji,jj,jk) ) / fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
 
@@ -382 +399 @@
                !
                DO jk = 1, nksr
-                  etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk)
+                  etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) )
                END DO
                etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero
@@ -404 +421 @@
                      zc0 = rn_abs * EXP( -fsdepw(ji,jj,jk  )*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk  )*zsi1r )
                      zc1 = rn_abs * EXP( -fsdepw(ji,jj,jk+1)*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk+1)*zsi1r )
-                     etot3(ji,jj,jk) = ro0cpr * (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) / fse3t(ji,jj,jk)
+                     etot3(ji,jj,jk) = ro0cpr * (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  )
                   END DO
                END DO