source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/limit_read_mod.F90

Last change on this file was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 18.0 KB
Line 
1!
2! $Header$
3!
4MODULE limit_read_mod
5!
6! This module reads the fichier "limit.nc" containing fields for surface forcing.
7!
8! Module subroutines :
9!  limit_read_frac    : call limit_read_tot and return the fractions
10!  limit_read_rug_alb : return rugosity and albedo, if coupled ocean call limit_read_tot first
11!  limit_read_sst     : return sea ice temperature   
12!  limit_read_tot     : read limit.nc and store the fields in local modules variables
13!
14  IMPLICIT NONE
15
16  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE, PRIVATE :: pctsrf
17!$OMP THREADPRIVATE(pctsrf)
18  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: rugos
19!$OMP THREADPRIVATE(rugos)
20  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: albedo
21!$OMP THREADPRIVATE(albedo) 
22  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: sst
23!$OMP THREADPRIVATE(sst) 
24#ifdef ISO
25  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: tuoce
26!$OMP THREADPRIVATE(tuoce)
27#endif
28  LOGICAL,SAVE :: read_continents=.FALSE.
29!$OMP THREADPRIVATE(read_continents)
30
31CONTAINS
32!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33!!
34!! Public subroutines :
35!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
36
37  SUBROUTINE limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified)
38!
39! This subroutine is called from "change_srf_frac" for case of
40! ocean=force or from ocean_slab_frac for ocean=slab.
41! The fraction for all sub-surfaces at actual time step is returned.
42
43    USE dimphy
44    USE indice_sol_mod
45
46! Input arguments
47!****************************************************************************************
48    INTEGER, INTENT(IN) :: itime   ! time step
49    INTEGER, INTENT(IN) :: jour    ! current day
50    REAL   , INTENT(IN) :: dtime   ! length of time step
51 
52! Output arguments
53!****************************************************************************************
54    REAL, DIMENSION(klon,nbsrf), INTENT(OUT) :: pctsrf_new  ! sub surface fractions
55    LOGICAL, INTENT(OUT)                     :: is_modified ! true if pctsrf is modified at this time step
56
57! End declaration
58!****************************************************************************************
59
60! 1) Read file limit.nc
61    CALL limit_read_tot(itime, dtime, jour, is_modified)
62
63! 2) Return the fraction read in limit_read_tot
64    pctsrf_new(:,:) = pctsrf(:,:)
65   
66  END SUBROUTINE limit_read_frac
67
68!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
69
70  SUBROUTINE limit_read_rug_alb(itime, dtime, jour, &
71       knon, knindex, &
72       rugos_out, alb_out &
73#ifdef ISO
74     & ,P_vegetation  &
75#endif     
76     & )
77!
78! This subroutine is called from surf_land_bucket.
79! The flag "ok_veget" must can not be true. If coupled run, "ocean=couple"
80! then this routine will call limit_read_tot.
81!
82    USE dimphy
83    USE surface_data
84#ifdef ISO
85    USE infotrac_phy, ONLY: niso
86    USE isotopes_mod, ONLY: P_veg
87#endif
88
89! Input arguments
90!****************************************************************************************
91    INTEGER, INTENT(IN) :: itime                     ! numero du pas de temps courant
92    INTEGER, INTENT(IN) :: jour                      ! jour a lire dans l'annee
93    REAL   , INTENT(IN) :: dtime                     ! pas de temps de la physique (en s)
94    INTEGER, INTENT(IN) :: knon                      ! nomber of points on compressed grid
95    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
96! Output arguments
97!****************************************************************************************
98    REAL, DIMENSION(klon), INTENT(OUT) :: rugos_out
99    REAL, DIMENSION(klon), INTENT(OUT) :: alb_out
100   
101! Local variables
102!****************************************************************************************
103    INTEGER :: i
104    LOGICAL :: is_modified
105
106
107#ifdef ISO 
108  real, intent(out), dimension(klon) :: P_vegetation
109#endif
110!****************************************************************************************
111
112IF (type_ocean == 'couple'.OR. &
113         (type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
114       ! limit.nc has not yet been read. Do it now!
115       CALL limit_read_tot(itime, dtime, jour, is_modified)
116    END IF
117
118    DO i=1,knon
119       rugos_out(i) = rugos(knindex(i))
120       alb_out(i)  = albedo(knindex(i))
121    END DO
122
123#ifdef ISO
124        !write(*,*) 'limit_read 120'
125        !* lecture de la fraction de l'évap sous forme de transpiration
126        !1) premier cas: on fixe P_vegetation à une valeur constante
127!        write(*,*) 'interfsurf 2990: P_veg=',P_veg
128        do i=1,knon
129           P_vegetation(i)=P_veg
130        enddo
131        ! 2) 2e cas: on lit une carte
132        ! A FAIRE
133
134        ! end verif
135#endif   
136
137  END SUBROUTINE limit_read_rug_alb
138
139!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140
141  SUBROUTINE limit_read_sst(knon, knindex, sst_out &
142#ifdef ISO
143     &  ,Roce,rlat   &
144#endif           
145    )
146!
147! This subroutine returns the sea surface temperature already read from limit.nc.
148! On rajoute la sortie tritium ocean sur le meme modele que la sst
149!
150    USE dimphy, ONLY : klon
151#ifdef ISO
152    USE infotrac_phy, ONLY: niso
153    USE isotopes_mod, ONLY: tcorr,toce,modif_sst, &
154   &    deltaTtest,sstlatcrit,deltaTtestpoles,dsstlatcrit, &
155   &    iso_HTO,ok_prod_nucl_tritium
156#ifdef ISOVERIF
157    USE isotopes_verif_mod, ONLY: iso_verif_egalite_vect2D,iso_verif_positif, &
158        iso_verif_positif_nostop
159#endif
160#endif
161
162    INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
163    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
164    REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out
165#ifdef ISO
166    REAL, DIMENSION(klon) :: tuoce_out ! sortie tritium surface ocean
167#endif
168
169    INTEGER :: i
170#ifdef ISO
171  real, intent(out), dimension(niso,klon) :: Roce
172  integer :: ixt 
173  real, intent(in),dimension(klon) :: rlat
174  real lat_locale
175#endif
176!#ifdef ISOVERIF
177!   integer iso_verif_positif_nostop
178!#endif
179
180    DO i = 1, knon
181       sst_out(i) = sst(knindex(i))
182    END DO
183
184#ifdef ISO
185     if (iso_HTO.gt.0) then
186     if (ok_prod_nucl_tritium) then ! si on active la production nucleaire de tritium
187        DO i = 1, knon
188          tuoce_out(i)=tuoce(knindex(i))
189        END DO
190     endif
191     endif
192#endif
193
194#ifdef ISO
195  if (modif_sst.ge.1) then
196  do i = 1, knon
197    lat_locale=rlat(knindex(i)) 
198    ! test: modification uniforme de la sst
199    if (modif_sst.eq.1) then   
200       sst_out(i)= sst_out(i)+deltaTtest 
201    elseif (modif_sst.eq.2) then   !if (modif_sst.eq.1) then       
202        ! pattern parabolique en dehors des tropiques (sstlatcrit)
203        if (abs(lat_locale).gt.sstlatcrit) then
204          sst_out(i)= sst_out(i)+deltaTtestpoles &
205   &             *(lat_locale**2-sstlatcrit**2) &
206   &             /(90.0**2-sstlatcrit**2)                     
207        endif !if (abs(lat_locale).gt.abs(sstlatcrit)) then
208
209    else if (modif_sst.eq.3) then
210
211        if (abs(lat_locale).gt.abs(sstlatcrit)) then
212            if (abs(lat_locale).gt.sstlatcrit+dsstlatcrit) then
213                sst_out(i)= sst_out(i)+deltaTtestpoles
214            else
215                sst_out(i)= sst_out(i)+deltaTtestpoles &
216    &               *(abs(lat_locale)-sstlatcrit)/dsstlatcrit
217            endif
218        endif
219    endif !if (modif_sst.eq.1) then   
220    enddo !do i = 1, knon
221    endif !if (modif_sst.ge.1) then
222#endif
223#ifdef ISOVERIF
224     do i=1,knon,20
225       call iso_verif_positif(sst_out(i)-100.0,'limit_read 4323')
226     enddo
227#endif
228
229
230#ifdef ISO
231        !* lecture de Roce
232        ! 1) première possibilité: valeur fixe à SMOW
233        DO i = 1, knon
234          do ixt=1,niso
235            Roce(ixt,i)=tcorr(ixt)*toce(ixt)
236          enddo !do ixt=1,niso
237        enddo !DO i = 1, knon
238        ! 2) deuxième possibilité: lecture de la carte
239        ! A FAIRE
240
241        ! lecture pour le tritium
242        if ((iso_HTO.gt.0).and.(ok_prod_nucl_tritium)) then
243             ! lecture de la carte tritium ocean surface
244             Roce(iso_HTO,i)=tcorr(iso_HTO)*tuoce_out(i)*1.E-18*2.
245        endif       
246#endif 
247
248#ifdef ISOVERIF
249        do i=1,knon
250          if (iso_verif_positif_nostop(370.0-sst_out(i), &
251              'limit_read 368').eq.1) then
252             write(*,*) 'i,knindex,sst_out=',i,knindex,sst_out(i)
253             stop
254          endif
255        enddo
256#endif     
257
258  END SUBROUTINE limit_read_sst
259
260!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
261!!
262!! Private subroutine :
263!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
264
265  SUBROUTINE limit_read_tot(itime, dtime, jour, is_modified)
266!
267! Read everything needed from limit.nc
268!
269! 0) Initialize
270! 1) Open the file limit.nc, if it is time
271! 2) Read fraction, if not type_ocean=couple
272! 3) Read sea surface temperature, if not type_ocean=couple
273! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
274! 5) Close file and distribuate variables to all processus
275
276    USE dimphy
277    USE mod_grid_phy_lmdz
278    USE mod_phys_lmdz_para
279    USE surface_data, ONLY : type_ocean, ok_veget
280    USE netcdf
281    USE indice_sol_mod
282#ifdef ISO
283    !USE infotrac_phy, ONLY: use_iso
284    USE isotopes_mod, ONLY : iso_HTO,ok_prod_nucl_tritium
285#ifdef ISOVERIF
286    USE isotopes_verif_mod, ONLY : iso_verif_positif_nostop
287#endif
288#endif
289
290    IMPLICIT NONE
291   
292! In- and ouput arguments
293!****************************************************************************************
294    INTEGER, INTENT(IN) :: itime   ! numero du pas de temps courant
295    INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
296    REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
297
298    LOGICAL, INTENT(OUT) :: is_modified  ! true if pctsrf is modified at this time step
299
300! Locals variables with attribute SAVE
301!****************************************************************************************
302! frequence de lecture des conditions limites (en pas de physique)
303    INTEGER,SAVE                              :: lmt_pas
304!$OMP THREADPRIVATE(lmt_pas)
305    LOGICAL, SAVE                             :: first_call=.TRUE.
306!$OMP THREADPRIVATE(first_call) 
307    INTEGER, SAVE                             :: jour_lu = -1
308!$OMP THREADPRIVATE(jour_lu) 
309! Locals variables
310!****************************************************************************************
311    INTEGER                                   :: nid, nvarid
312    INTEGER                                   :: ii, ierr
313    INTEGER, DIMENSION(2)                     :: start, epais
314    REAL, DIMENSION(klon_glo,nbsrf)           :: pct_glo  ! fraction at global grid
315    REAL, DIMENSION(klon_glo)                 :: sst_glo  ! sea-surface temperature at global grid
316    REAL, DIMENSION(klon_glo)                 :: rug_glo  ! rugosity at global grid
317    REAL, DIMENSION(klon_glo)                 :: alb_glo  ! albedo at global grid
318    CHARACTER(len=20)                         :: modname='limit_read_mod'
319#ifdef ISO
320    REAL, DIMENSION(klon_glo)                 :: tuoce_glo  ! sea-surface tritium et global grid
321#endif     
322
323! End declaration
324!****************************************************************************************
325
326!****************************************************************************************
327! 0) Initialization
328!
329!****************************************************************************************
330    IF (first_call) THEN
331       ! calculate number of time steps for one day
332       lmt_pas = NINT(86400./dtime * 1.0)
333       
334       ! Allocate module save variables
335       IF ( type_ocean /= 'couple' ) THEN
336          ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
337          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1)
338       END IF
339
340       IF ( .NOT. ok_veget ) THEN
341          ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
342          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating rugos and albedo',1)
343       END IF
344
345       first_call=.FALSE.
346    ENDIF
347 
348!****************************************************************************************
349! 1) Open the file limit.nc if it is the right moment to read, once a day.
350!    The file is read only by the master thread of the master mpi process(is_mpi_root)
351!
352!****************************************************************************************
353
354    is_modified = .FALSE.
355    IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN   ! time to read
356       jour_lu = jour
357       is_modified = .TRUE.
358!$OMP MASTER  ! Only master thread
359       IF (is_mpi_root) THEN ! Only master processus
360
361          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
362          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
363               'Pb d''ouverture du fichier de conditions aux limites',1)
364         
365          ! La tranche de donnees a lire:
366          start(1) = 1
367          start(2) = jour
368          epais(1) = klon_glo
369          epais(2) = 1
370
371
372!****************************************************************************************
373! 2) Read fraction if not type_ocean=couple
374!
375!****************************************************************************************
376
377          IF ( type_ocean /= 'couple') THEN
378!
379! Ocean fraction
380             ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
381             IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
382             
383             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
384             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
385!
386! Sea-ice fraction
387             ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
388             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
389
390             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
391             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
392
393
394! Read land and continentals fraction only if asked for
395             IF (read_continents .OR. itime == 1) THEN
396!
397! Land fraction
398                ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
399                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
400               
401                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
402                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
403!
404! Continentale ice fraction
405                ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
406                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
407
408                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
409                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
410             END IF
411
412          END IF ! type_ocean /= couple
413
414!****************************************************************************************
415! 3) Read sea-surface temperature, if not coupled ocean
416!
417!****************************************************************************************
418          IF ( type_ocean /= 'couple') THEN
419
420             ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
421             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
422
423             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
424             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
425
426#ifdef ISO
427             IF ((iso_HTO.gt.0).and.(ok_prod_nucl_tritium)) THEN
428               ierr = NF90_INQ_VARID(nid, 'TUOCE', nvarid)
429               IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Le champ <TUOCE> est absent',1)
430
431               ierr = NF90_GET_VAR(nid,nvarid,tuoce_glo,start,epais)
432               IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <TUOCE>',1)
433             END IF
434#ifdef ISOVERIF
435             do ii=1,klon_glo
436               if (iso_verif_positif_nostop(370.0-sst_glo(ii),  &
437                'limit_read 384').eq.1) then
438                 write(*,*) 'ii,sst_glo=',ii,sst_glo(ii)
439                 write(*,*) 'jour,start,epais=',jour,start,epais
440                 stop
441               endif
442             enddo
443#endif
444#endif
445
446          END IF
447
448!****************************************************************************************
449! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
450!
451!****************************************************************************************
452
453          IF (.NOT. ok_veget) THEN
454!
455! Read albedo
456             ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
457             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
458
459             ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
460             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
461!
462! Read rugosity
463             ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
464             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
465
466             ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
467             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
468
469          END IF
470
471!****************************************************************************************
472! 5) Close file and distribuate variables to all processus
473!
474!****************************************************************************************
475          ierr = NF90_CLOSE(nid)
476          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
477       ENDIF ! is_mpi_root
478
479!$OMP END MASTER
480!$OMP BARRIER
481
482       IF ( type_ocean /= 'couple') THEN
483          CALL Scatter(sst_glo,sst)
484          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
485          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
486          IF (read_continents .OR. itime == 1) THEN
487             CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
488             CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
489          END IF
490#ifdef ISO
491          IF ((iso_HTO.gt.0).and.(ok_prod_nucl_tritium)) THEN
492             CALL Scatter(tuoce_glo,tuoce)
493          END IF
494#endif
495       END IF
496
497       IF (.NOT. ok_veget) THEN
498          CALL Scatter(alb_glo, albedo)
499          CALL Scatter(rug_glo, rugos)
500       END IF
501
502    ENDIF ! time to read
503
504  END SUBROUTINE limit_read_tot
505
506
507END MODULE limit_read_mod
Note: See TracBrowser for help on using the repository browser.