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.6_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90 @ 7803

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

Merge in changes to allow a maximum chlorophyll assimilation increment to be specified. See internal UKMO ticket 690.
File size: 11.8 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
29CONTAINS
30
31   SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
32      &                             k_maxchlinc, logchl_bkginc, logchl_balinc )
33      !!---------------------------------------------------------------------------
34      !!                    ***  ROUTINE asm_logchl_bal_ersem  ***
35      !!
36      !! ** Purpose :   calculate increments to ERSEM from logchl increments
37      !!
38      !! ** Method  :   convert logchl increments to chl increments
39      !!                split between the ERSEM PFTs
40      !!                spread through the mixed layer
41      !!                [forthcoming: calculate increments to nutrients and zooplankton]
42      !!
43      !! ** Action  :   populate logchl_balinc
44      !!
45      !! References :   forthcoming...
46      !!---------------------------------------------------------------------------
47      !!
48      LOGICAL,  INTENT(in   )                               :: ld_logchlpftinc
49      INTEGER,  INTENT(in   )                               :: npfts
50      INTEGER,  INTENT(in   )                               :: mld_choice_bgc
51      REAL(wp), INTENT(in   )                               :: k_maxchlinc
52      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts)     :: logchl_bkginc
53      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc
54      !!
55      INTEGER                      :: ji, jj, jk
56      INTEGER                      :: jkmax
57      REAL(wp)                     :: chl_tot, chl_inc
58      REAL(wp), DIMENSION(jpi,jpj) :: zmld
59      !!---------------------------------------------------------------------------
60
61      ! Split surface logchl incs into surface Chl1-4 incs
62      !
63      ! In order to transform logchl incs to chl incs, need to account for the background,
64      ! cannot simply do 10^logchl_bkginc. Need to:
65      ! 1) Add logchl inc to log10(background) to get log10(analysis)
66      ! 2) Take 10^log10(analysis) to get analysis
67      ! 3) Subtract background from analysis to get chl incs
68      ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
69      !
70      ! Only apply increments if all of Chl1-4 background values are > 0
71      ! In theory, they always will be, and if any are not that's a sign
72      ! that something's going wrong which the assimilation might make worse
73      !
74      IF ( ld_logchlpftinc ) THEN
75         !
76         ! Assimilating separate PFTs, so separately transform each from LogChl to Chl
77         !
78         IF ( npfts /= 4 ) THEN
79            CALL ctl_stop( 'If assimilating PFTs into ERSEM, nn_asmpfts must be 4' )
80         ENDIF
81         DO jj = 1, jpj
82            DO ji = 1, jpi
83               IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
84                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
85                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
86                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
87                  IF ( logchl_bkginc(ji,jj,1) /= 0.0 ) THEN
88                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) ) + &
89                        &                                                   logchl_bkginc(ji,jj,1) )                      - &
90                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
91                  ENDIF
92                  IF ( logchl_bkginc(ji,jj,2) /= 0.0 ) THEN
93                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) ) + &
94                        &                                                   logchl_bkginc(ji,jj,2) )                      - &
95                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
96                  ENDIF
97                  IF ( logchl_bkginc(ji,jj,3) /= 0.0 ) THEN
98                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) ) + &
99                        &                                                   logchl_bkginc(ji,jj,3) )                      - &
100                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
101                  ENDIF
102                  IF ( logchl_bkginc(ji,jj,4) /= 0.0 ) THEN
103                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) + &
104                        &                                                   logchl_bkginc(ji,jj,4) )                      - &
105                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
106                  ENDIF
107                  IF (k_maxchlinc > 0.0) THEN
108                     chl_inc = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
109                               logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
110                               logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
111                               logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
112                     IF ( ABS(chl_inc) > k_maxchlinc ) THEN
113                        chl_tot = ABS(chl_inc) / k_maxchlinc
114                        logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / chl_tot
115                        logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / chl_tot
116                        logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / chl_tot
117                        logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / chl_tot
118                     ENDIF
119                  ENDIF
120               ENDIF
121            END DO
122         END DO
123      ELSE
124         !
125         ! Assimilating total Chl, so transform total from LogChl to Chl
126         ! and split between PFTs according to the existing background ratios
127         !
128         IF ( npfts /= 1 ) THEN
129            CALL ctl_stop( 'If assimilating total chlorophyll, nn_asmpfts must be 1' )
130         ENDIF
131         DO jj = 1, jpj
132            DO ji = 1, jpi
133               IF ( ( logchl_bkginc(ji,jj,1)               /= 0.0 ) .AND. &
134                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) >  0.0 ) .AND. &
135                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) >  0.0 ) .AND. &
136                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) >  0.0 ) .AND. &
137                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) >  0.0 ) ) THEN
138                  chl_tot = trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
139                     &      trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
140                  chl_inc = 10**( LOG10( chl_tot ) + logchl_bkginc(ji,jj,1) ) - chl_tot
141                  IF (k_maxchlinc > 0.0) THEN
142                     chl_inc = MAX( -1.0 * k_maxchlinc, MIN( chl_inc, k_maxchlinc ) )
143                  ENDIF
144                  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
145                  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
146                  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
147                  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
148               ENDIF
149            END DO
150         END DO
151      ENDIF
152
153      ! Propagate surface Chl1-4 incs through mixed layer
154      ! First, choose mixed layer definition
155      !
156      SELECT CASE( mld_choice_bgc )
157      CASE ( 1 )                   ! Turbocline/mixing depth [W points]
158         zmld(:,:) = hmld(:,:)
159      CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
160         zmld(:,:) = hmlp(:,:)
161      CASE ( 3 )                   ! Kara MLD [Interpolated]
162#if defined key_karaml
163         IF ( ln_kara ) THEN
164            zmld(:,:) = hmld_kara(:,:)
165         ELSE
166            CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
167               &           ' but ln_kara=.false.' )
168         ENDIF
169#else
170         CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
171            &           ' but is not defined' )
172#endif
173      CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
174         zmld(:,:) = hmld_tref(:,:)
175      CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
176         zmld(:,:) = hmlpt(:,:)
177      END SELECT
178      !
179      ! Now set MLD to bottom of a level and propagate incs equally through mixed layer
180      !
181      DO jj = 1, jpj
182         DO ji = 1, jpi
183            !
184            jkmax = jpk-1
185            DO jk = jpk-1, 1, -1
186               IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
187                  & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
188                  zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
189                  jkmax = jk
190               ENDIF
191            END DO
192            !
193            DO jk = 2, jkmax
194               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
195               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
196               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
197               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
198            END DO
199            !
200         END DO
201      END DO
202
203      ! Multivariate balancing forthcoming...
204
205   END SUBROUTINE asm_logchl_bal_ersem
206
207#else
208   !!----------------------------------------------------------------------
209   !!   Default option : Empty routine
210   !!----------------------------------------------------------------------
211CONTAINS
212   SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
213      &                             k_maxchlinc, logchl_bkginc, logchl_balinc )
214      LOGICAL :: ld_logchlpftinc
215      INTEGER :: npfts
216      INTEGER :: mld_choice_bgc
217      REAL    :: k_maxchlinc
218      REAL    :: logchl_bkginc(:,:,:)
219      REAL    :: logchl_balinc(:,:,:,:)
220      WRITE(*,*) 'asm_logchl_bal_ersem: You should not have seen this print! error?'
221   END SUBROUTINE asm_logchl_bal_ersem
222#endif
223
224   !!======================================================================
225END MODULE asmlogchlbal_ersem
Note: See TracBrowser for help on using the repository browser.