source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 2437

Last change on this file since 2437 was 2311, checked in by emillour, 6 years ago

Mars GCM:
Code tidying: use getin_p() instead of getin() and use "call abort_physic"
instead of "stop" or "call abort".
EM
File size: 17.1 KB
Line 
1      SUBROUTINE surfini(ngrid,piceco2,qsurf)
2
3      USE ioipsl_getin_p_mod, ONLY : getin_p
4      use netcdf
5      use tracer_mod, only: nqmx, noms
6      use geometry_mod, only: longitude, latitude ! in radians
7      use surfdat_h, only: watercaptag, frost_albedo_threshold,
8     &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
9     &                     albedice, dryness
10#ifndef MESOSCALE
11      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
12      use mod_phys_lmdz_para, only : is_master, gather, scatter
13#endif
14      USE comcstfi_h, ONLY: pi
15      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
16      use datafile_mod, only: datadir
17      IMPLICIT NONE
18c=======================================================================
19c
20c   creation des calottes pour l'etat initial
21c
22c=======================================================================
23c-----------------------------------------------------------------------
24c   Declarations:
25c   -------------
26      include "callkeys.h"
27
28      integer,intent(in) :: ngrid ! number of atmospheric columns
29      real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
30      real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
31
32      INTEGER ig,icap,iq,alternate
33      REAL icedryness ! ice dryness
34     
35      ! longwatercaptag is watercaptag. Trick for some compilers
36      LOGICAL, DIMENSION(100000) :: longwatercaptag
37     
38! There are 3 different modes for ice distribution:
39! icelocationmode = 1 ---> based on data from surface.nc
40! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
41! icelocationmode = 3 ---> based on logical relations for latitude and longitude
42! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
43      INTEGER,SAVE :: icelocationmode = 2
44       
45       
46      !in case icelocationmode == 1
47      INTEGER i,j
48      INTEGER     imd,jmd
49      PARAMETER   (imd=360,jmd=180)
50      REAL        zdata(imd,jmd)
51      REAL        zelat,zelon
52
53#ifndef MESOSCALE
54      INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
55#endif
56      INTEGER latice(nbp_lat-1,2),lonice (nbp_lon,2) ! number of counts | detected ice along lat & lon axis
57
58      REAL step,count,ratiolat
59
60      INTEGER   ierr,nid,nvarid
61     
62      REAL,SAVE :: min_icevalue = 500.
63      character(len=50) :: string = 'thermal'
64     
65      character (len=100) :: zedatafile
66
67#ifdef MESOSCALE
68
69      do ig=1,ngrid
70
71         !write(*,*) "all qsurf to zero. dirty."
72         do iq=1,nqmx
73         qsurf(ig,iq)=0.  !! on jette les inputs GCM
74                          !! on regle juste watercaptag
75                          !! il faudrait garder les inputs GCM
76                          !! si elles sont consequentes
77         enddo
78         if ( ( latitude(ig)*180./pi .gt. 70. ) .and.
79     .        ( albedodat(ig) .ge. 0.26   ) )  then
80                 write(*,*)"outlier ",ig
81                 watercaptag(ig)  = .true.
82                 dryness(ig)      = 1.
83                 albedodat(ig)    = albedo_h2o_ice  !! pour output
84         else
85                 watercaptag(ig)  = .false.
86                 dryness(ig)      = 1.
87         endif
88
89      enddo
90#endif
91! problem with nested precompiling flags
92
93#ifndef MESOSCALE
94      ! to handle parallel cases
95#if CPP_PARA
96      logical watercaptag_glo(klon_glo)
97      real dryness_glo(klon_glo)
98      real lati_glo(klon_glo)
99      real long_glo(klon_glo)
100#else
101      logical watercaptag_glo(ngrid)
102      real dryness_glo(ngrid)
103      real lati_glo(ngrid)
104      real long_glo(ngrid)
105#endif
106#endif
107
108#ifndef MESOSCALE
109
110c
111c=======================================================================
112! Initialize watercaptag (default is false)
113      watercaptag_glo(:)=.false.
114
115c     water ice outliers
116c     ------------------------------------------
117
118      IF ((water) .and. (caps)) THEN
119     
120c Perennial H20 north cap defined by watercaptag=true (allows surface to be
121c hollowed by sublimation in vdifc).
122
123c We might not want albedodat to be modified because it is used to write
124c restart files. Instead, albedo is directly modified when needed (i.e.
125c if we have watercaptag and no co2 ice), below and in albedocaps.F90
126
127c       "Dryness coefficient" controlling the evaporation and
128c        sublimation from the ground water ice (close to 1)
129c        HERE, the goal is to correct for the fact
130c        that the simulated permanent water ice polar caps
131c        is larger than the actual cap and the atmospheric
132c        opacity not always realistic.
133
134         alternate = 0
135         
136         if (ngrid .ne. 1) then
137           watercaptag(:) = .false.
138           longwatercaptag(:) = .false.
139         endif
140         
141         write(*,*) "surfini: Ice dryness ?"
142         icedryness=1. ! default value
143         call getin_p("icedryness",icedryness)
144         write(*,*) "surfini: icedryness = ",icedryness
145         dryness (:) = icedryness
146         
147      ! To be able to run in parallel, we work on the full grid
148      ! and dispatch results afterwards
149
150      ! start by geting latitudes and logitudes on full grid
151      ! (in serial mode, this is just a copy)
152      call gather(latitude,lati_glo)
153      call gather(longitude,long_glo)
154
155      if (is_master) then
156
157        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
158     
159         print*, 'ngrid = 1, do no put ice caps in surfini.F'
160
161        ELSE IF (icelocationmode .eq. 1) THEN
162     
163         print*,'Surfini: ice caps defined from surface.nc'
164           
165! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
166! Typically, it is for thermal inertia above 500 tiu.
167! Two conditions are verified:
168! 1. GCM ice caps are defined such as area is conserved for a given latitude
169! (the approximation is that all points within the GCM latitude resolution have the same area).
170! 2. caps are placed to fill the GCM points with the most detected ice first.
171     
172
173           
174         zedatafile = trim(datadir)
175 
176       
177         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
178     &   NF90_NOWRITE,nid)
179     
180         IF (ierr.NE.nf90_noerr) THEN
181       write(*,*)'Error : cannot open file surface.nc '
182       write(*,*)'(in phymars/surfini.F)'
183       write(*,*)'It should be in :',trim(zedatafile),'/'
184       write(*,*)'1) You can set this path in the callphys.def file:'
185       write(*,*)'   datadir=/path/to/the/datafiles'
186       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
187       write(*,*)'   can be obtained online on:'
188       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
189       call abort_physic("surfini","missing surface.nc file",1)
190         ENDIF
191     
192     
193         ierr=nf90_inq_varid(nid, string, nvarid)
194         if (ierr.ne.nf90_noerr) then
195          write(*,*) 'surfini error, cannot find ',trim(string)
196          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
197          write(*,*)trim(nf90_strerror(ierr))
198          call abort_physic("surfini","missing "//trim(string),1)
199         endif
200
201         ierr=nf90_get_var(nid, nvarid, zdata)
202
203         if (ierr.ne.nf90_noerr) then
204          write(*,*) 'surfini: error failed loading ',trim(string)
205          write(*,*)trim(nf90_strerror(ierr))
206          call abort_physic("surfini","failed loading "//trim(string),1)
207         endif
208 
209                     
210         ierr=nf90_close(nid)
211 
212
213         nb_ice(:,1) = 1 ! default: there is no ice
214         latice(:,1) = 1
215         lonice(:,1) = 1
216         nb_ice(:,2) = 0
217         latice(:,2) = 0
218         lonice(:,2) = 0
219         !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
220
221         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
222         do ig=2,klon_glo-1
223     
224          ! loop over the surface file grid     
225          do i=1,imd
226           do j=1,jmd
227             zelon = i - 180.
228             zelat = 90. - j
229            if ((abs(lati_glo(ig)*180./pi-zelat).le.
230     &           90./real(nbp_lat-1)) .and.
231     &          (abs(long_glo(ig)*180./pi-zelon).le.
232     &           180./real(nbp_lon))) then
233              ! count all points in that GCM grid point
234              nb_ice(ig,1) = nb_ice(ig,1) + 1
235              if (zdata(i,j) > min_icevalue)
236                 ! count all detected points in that GCM grid point
237     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
238             endif
239           enddo
240          enddo 
241
242        ! projection of nb_ice on GCM lat and lon axes
243          latice(1+(ig-2)/nbp_lon,:) =
244     &     latice(1+(ig-2)/nbp_lon,:) + nb_ice(ig,:)
245          lonice(1+mod(ig-2,nbp_lon),:) =
246     &     lonice(1+mod(ig-2,nbp_lon),:) + nb_ice(ig,:) ! lonice is USELESS ...
247
248         enddo ! of do ig=2,klon_glo-1
249     
250
251     
252         ! special case for poles
253         nb_ice(1,2)   = 1  ! ice prescribed on north pole
254         latice(1,:)   = nb_ice(1,:)
255         lonice(1,:)   = nb_ice(1,:)
256         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
257         lonice(nbp_lon,:) = nb_ice(ngrid,:)
258     
259     
260!      print*, 'latice TOT', latice(:,1)
261!      print*, 'latice FOUND', latice(:,2)
262!      print*, 'lonice TOT', lonice(:,1)
263!      print*, 'lonice FOUND', lonice(:,2)
264     
265!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
266!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
267     
268!      print*,''
269!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
270!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
271     
272   
273         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
274         do i=1,(nbp_lat-1)/2
275          step  = 1. ! threshold to add ice cap
276          count = 0. ! number of ice GCM caps at this latitude
277          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
278          ratiolat  = real(latice(i,2))/real(latice(i,1))
279          !print*,'i',i,(i-1)*iim+2,i*iim+1
280     
281          ! put ice caps while there is not enough ice,
282          ! as long as the threshold is above 20%
283          do while ((count.le.ratiolat*nbp_lon).and.(step.ge.0.2))
284           count = 0.
285           ! loop over GCM longitudes
286           do j=1,nbp_lon
287            ! if the detected ice ratio in the GCM grid point
288            ! is more than 'step', then add ice
289            if (real(nb_ice((i-1)*nbp_lon+1+j,2))
290     &        / real(nb_ice((i-1)*nbp_lon+1+j,1)) .ge. step) then
291                  watercaptag_glo((i-1)*nbp_lon+1+j) = .true.
292                  count = count + 1
293            endif
294           enddo ! of do j=1,nbp_lon
295           !print*, 'step',step,count,ratiolat*nbp_lon
296           step = step - 0.01
297          enddo ! of do while
298          !print*, 'step',step,count,ratiolat*nbp_lon
299
300         enddo ! of do i=1,jjm/2
301           
302
303        ELSE IF (icelocationmode .eq. 2) THEN
304     
305         print*,'Surfini: predefined ice caps'
306     
307         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
308           
309          print*,'water ice caps distribution for 32x24 resolution'
310          longwatercaptag(1:9)    = .true. ! central cap - core
311          longwatercaptag(26:33)  = .true. ! central cap
312          longwatercaptag(1:33)  = .true. ! central cap
313          longwatercaptag(56)  = .true. ! central cap
314          longwatercaptag(58)  = .true. ! central cap
315          longwatercaptag(60)  = .true. ! central cap
316          longwatercaptag(62)  = .true. ! central cap
317          longwatercaptag(64)  = .true. ! central cap
318!---------------------   OUTLIERS  ----------------------------
319
320         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48
321
322          print*,'water ice caps distribution for 64x48 resolution'
323          longwatercaptag(1:65)   = .true. ! central cap - core
324          longwatercaptag(75:85)  = .true. ! central cap
325          longwatercaptag(93:114) = .true. ! central cap
326!---------------------   OUTLIERS  ----------------------------
327          if (.true.) then
328          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
329          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
330          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
331          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
332          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
333          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
334          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
335          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
336          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
337          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
338          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
339          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
340          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
341          longwatercaptag(194)    = .true. ! outlier, lat = 75
342          longwatercaptag(203)    = .true. ! outlier, lat = 75
343          longwatercaptag(207)    = .true. ! outlier, lat = 75
344          longwatercaptag(244)    = .true. ! outlier, lat = 75
345          longwatercaptag(246)    = .true. ! outlier, lat = 75
346          longwatercaptag(250)    = .true. ! outlier, lat = 75
347          longwatercaptag(252)    = .true. ! outlier, lat = 75
348          longwatercaptag(254)    = .true. ! outlier, lat = 75
349          longwatercaptag(256)    = .true. ! outlier, lat = 75
350          endif
351!--------------------------------------------------------------       
352
353
354           
355         else if (klon_glo .ne. 1) then
356       
357          print*,'No predefined ice location for this resolution :',
358     &           nbp_lon,nbp_lat-1
359          print*,'Please change icelocationmode in surfini.F'
360          print*,'Or add some new definitions ...'
361          call abort_physic("surfini",
362     &         "no pre-definitions for this resolution",1)
363         
364         endif
365
366         do ig=1,klon_glo
367          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
368         enddo
369
370
371        ELSE IF (icelocationmode .eq. 3) THEN
372     
373         print*,'Surfini: ice caps defined by lat and lon values'
374
375         do ig=1,klon_glo
376         
377c-------- Towards olympia planitia water caps -----------
378c--------------------------------------------------------
379
380          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
381     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
382     .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
383     .           ( long_glo(ig)*180./pi .le. 181. ) )
384     .         .or.
385
386     .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
387     .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
388     .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
389     .           ( long_glo(ig)*180./pi .le. 168. ) )
390     .         .or.
391     .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
392     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
393     .           ( long_glo(ig)*180./pi .ge. -150.) .and.
394     .           ( long_glo(ig)*180./pi .le. -110.) ) )
395     .         then
396             
397               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
398              !    watercaptag(ig)=.true.
399                  alternate = 1
400               else
401                  alternate = 0
402               endif !end if alternate = 0
403               
404          endif
405
406c----------- Opposite olympia planitia water cap --------
407c--------------------------------------------------------
408
409          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
410     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
411     .         .and.
412     .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
413     .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
414!     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
415!     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
416!     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
417!     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
418        !   watercaptag_glo(ig)=.true.
419          endif
420
421
422c -------------------- Central cap ----------------------
423c--------------------------------------------------------
424
425          if (abs(lati_glo(ig)*180./pi).gt.80)
426     .          watercaptag_glo(ig)=.true.
427           
428c--------------------------------------------------------
429c--------------------------------------------------------
430         end do ! of (klon_glo)
431
432
433        ELSE
434     
435         print*, 'In surfini.F, icelocationmode is ', icelocationmode
436         print*, 'It should be 1, 2 or 3.'
437         call abort_physic("surfini","wrong icelocationmode",1)
438
439        ENDIF ! of if (icelocation)
440       
441       
442        ! print caps locations - useful for plots too
443        print*,'surfini: latitude | longitude | ig'
444        do ig=1,klon_glo
445          dryness_glo(ig) = icedryness
446
447          if (watercaptag_glo(ig)) then
448             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
449     &              long_glo(ig)*180./pi, ig
450          endif
451        enddo
452       
453       endif !of if (is_master)
454       
455       if (ngrid.gt.1) then
456        ! Now scatter fields watercaptag and dryness from master to all
457        ! (is just a plain copy in serial mode)
458        call scatter(dryness_glo,dryness)
459        call scatter(watercaptag_glo,watercaptag)
460       endif
461       ELSE
462         watercaptag(:) = .false.
463       ENDIF ! (caps & water)
464! end of #else of #ifndef MESOSCALE
465#endif       
466
467      END
Note: See TracBrowser for help on using the repository browser.