New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asmlogchlbal_ersem.F90 in branches/UKMO/dev_r5518_v3.4_asm_nemovar_community_bgc_ersem/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_v3.4_asm_nemovar_community_bgc_ersem/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90 @ 7716

Last change on this file since 7716 was 7716, checked in by dford, 7 years ago

Tidy up and simplify use of arrays for logchl assimilation balancing.
File size: 9.6 KB
Line 
1MODULE asmlogchlbal_ersem
2   !!======================================================================
3   !!                       ***  MODULE asmlogchlbal_ersem  ***
4   !! Calculate increments to ERSEM based on surface logchl increments
5   !!======================================================================
6   !! History :  3.6  ! 2016-09  (D. Ford)     Original code
7   !!----------------------------------------------------------------------
8#if defined key_asminc && defined key_fabm
9   !!----------------------------------------------------------------------
10   !!   'key_asminc'         : assimilation increment interface
11   !!   'key_fabm'           : FABM-ERSEM coupling
12   !!----------------------------------------------------------------------
13   !!   asm_logchl_bal_ersem : routine to calculate increments to ERSEM
14   !!----------------------------------------------------------------------
15   USE par_kind,    ONLY: wp             ! kind parameters
16   USE par_oce,     ONLY: jpi, jpj, jpk  ! domain array sizes
17   USE dom_oce,     ONLY: gdepw_n        ! domain information
18   USE zdfmxl                            ! mixed layer depth
19   USE iom                               ! i/o
20   USE trc,         ONLY: trn            ! ERSEM state variables
21   USE par_fabm                          ! FABM parameters
22   USE par_trc,     ONLY: jptra          ! Tracer parameters
23
24   IMPLICIT NONE
25   PRIVATE                   
26
27   PUBLIC asm_logchl_bal_ersem
28
29   REAL(wp), PARAMETER :: tiny = 1e-20
30
31CONTAINS
32
33   SUBROUTINE asm_logchl_bal_ersem( ln_logchlpftinc, npfts, mld_choice_bgc, &
34      &                             logchl_bkginc, logchl_balinc )
35      !!---------------------------------------------------------------------------
36      !!                    ***  ROUTINE asm_logchl_bal_ersem  ***
37      !!
38      !! ** Purpose :   calculate increments to ERSEM from logchl increments
39      !!
40      !! ** Method  :   forthcoming...
41      !!
42      !! ** Action  :   forthcoming...
43      !!
44      !! References :   forthcoming...
45      !!---------------------------------------------------------------------------
46      !!
47      LOGICAL,  INTENT(in   )                               :: ln_logchlpftinc
48      INTEGER,  INTENT(in   )                               :: npfts
49      INTEGER,  INTENT(in   )                               :: mld_choice_bgc
50      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts)     :: logchl_bkginc
51      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc
52      !!
53      INTEGER                      :: ji, jj, jk
54      INTEGER                      :: jkmax
55      REAL(wp)                     :: chl_tot, chl_inc
56      REAL(wp), DIMENSION(jpi,jpj) :: zmld
57      !!---------------------------------------------------------------------------
58
59      ! Split surface logchl incs into surface Chl1-4 incs
60      !
61      ! In order to transform logchl incs to chl incs, need to account for the background,
62      ! cannot simply do 10^logchl_bkginc. Need to:
63      ! 1) Add logchl inc to log10(background) to get log10(analysis)
64      ! 2) Take 10^log10(analysis) to get analysis
65      ! 3) Subtract background from analysis to get chl incs
66      !
67      IF ( ln_logchlpftinc ) THEN
68         !
69         ! Assimilating separate PFTs, so separately transform each from LogChl to Chl
70         !
71         IF ( npfts /= 4 ) THEN
72            CALL ctl_stop( 'If assimilating PFTs into ERSEM, nn_asmpfts must be 4' )
73         ENDIF
74         DO jj = 1, jpj
75            DO ji = 1, jpi
76               IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
77                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
78                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
79                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
80                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) ) + &
81                     &                                                   logchl_bkginc(ji,jj,1) )                      - &
82                     &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
83                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) ) + &
84                     &                                                   logchl_bkginc(ji,jj,2) )                      - &
85                     &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
86                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) ) + &
87                     &                                                   logchl_bkginc(ji,jj,3) )                      - &
88                     &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
89                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) + &
90                     &                                                   logchl_bkginc(ji,jj,4) )                      - &
91                     &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
92               ENDIF
93            END DO
94         END DO
95      ELSE
96         !
97         ! Assimilating total Chl, so transform total from LogChl to Chl
98         ! and split between PFTs according to the existing background ratios
99         !
100         IF ( npfts /= 1 ) THEN
101            CALL ctl_stop( 'If assimilating total chlorophyll, nn_asmpfts must be 1' )
102         ENDIF
103         DO jj = 1, jpj
104            DO ji = 1, jpi
105               IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
106                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
107                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
108                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
109                  chl_tot = trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
110                     &      trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
111                  chl_inc = 10**( LOG10( chl_tot ) + logchl_bkginc(ji,jj,1) ) - chl_tot
112                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / chl_tot
113                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / chl_tot
114                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / chl_tot
115                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / chl_tot
116               ENDIF
117            END DO
118         END DO
119      ENDIF
120
121      ! Propagate surface Chl1-4 incs through mixed layer
122      ! First, choose mixed layer definition
123      !
124      SELECT CASE( mld_choice_bgc )
125      CASE ( 1 )                   ! Turbocline/mixing depth [W points]
126         zmld(:,:) = hmld(:,:)
127      CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
128         zmld(:,:) = hmlp(:,:)
129      CASE ( 3 )                   ! Kara MLD [Interpolated]
130#if defined key_karaml
131         IF ( ln_kara ) THEN
132            zmld(:,:) = hmld_kara(:,:)
133         ELSE
134            CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
135               &           ' but ln_kara=.false.' )
136         ENDIF
137#else
138         CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
139            &           ' but is not defined' )
140#endif
141      CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
142         zmld(:,:) = hmld_tref(:,:)
143      CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
144         zmld(:,:) = hmlpt(:,:)
145      END SELECT
146      !
147      ! Now set MLD to bottom of a level and propagate incs equally through mixed layer
148      !
149      DO jj = 1, jpj
150         DO ji = 1, jpi
151            !
152            jkmax = jpk-1
153            DO jk = jpk-1, 1, -1
154               IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
155                  & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
156                  zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
157                  jkmax = jk
158               ENDIF
159            END DO
160            !
161            DO jk = 2, jkmax
162               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
163               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
164               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
165               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
166            END DO
167            !
168         END DO
169      END DO
170
171      ! Multivariate balancing forthcoming...
172
173   END SUBROUTINE asm_logchl_bal_ersem
174
175#else
176   !!----------------------------------------------------------------------
177   !!   Default option : Empty routine
178   !!----------------------------------------------------------------------
179CONTAINS
180   SUBROUTINE asm_logchl_bal_ersem( ln_logchlpftinc, npfts, mld_choice_bgc, &
181      &                             logchl_bkginc, logchl_balinc )
182      LOGICAL :: ln_logchlpftinc
183      INTEGER :: npfts
184      INTEGER :: mld_choice_bgc
185      REAL    :: logchl_bkginc(:,:,:)
186      REAL    :: logchl_balinc(:,:,:,:)
187      WRITE(*,*) 'asm_logchl_bal_ersem: You should not have seen this print! error?'
188   END SUBROUTINE asm_logchl_bal_ersem
189#endif
190
191   !!======================================================================
192END MODULE asmlogchlbal_ersem
Note: See TracBrowser for help on using the repository browser.